Postrun data process script - 1PFT, 2sites, CO2 treatment case - CO2 treatment x site; instances - different hydraulic traits/parameters 1) import modules 2) pre-analysis a) read flux tower data b) read in lookup table c) read in selected history fields d) define functions to plot for stoma, xylem, and root profile 3) run data analysis a) total basal area, DDBH, mortality as BA(MTba), recruitment, GPP for each PFT (pine) b) calculate mean and std of BA, DDBH, MTba for each case and instance c) calculate difference in mean and std for each instance between different CO2 treatment cases 4) make plots a) choose instances to compare b) plot stoma, xylem, root profiles for selected instances c) plot
======Available history fields======================= hist_fincl1='NEP','NPP','GPP','QVEGT','EFLX_LH_TOT','AR','HR','ED_biomass','FSDS','FSH','Rnet','BTRAN','TLAI','ED_NPATCHES','AREA_PLANT','AREA_TREES','CANOPY_SPREAD','PFTbiomass','PFTcrownarea','RECRUITMENT','NCL_BY_AGE','QVEGE','ED_balive','DDBH_SCPF','BA_SCPF','NPLANT_SCPF','NPP_SCPF', 'CANOPY_HEIGHT_DIST','LEAF_HEIGHT_DIST','C_STOMATA_BY_AGE','PARPROF_DIR_CNLFPFT','CROWNAREA_CNLF','CROWNAREA_CAN','BA_SCLS','GPP_SCPF','ED_bsapwood','ED_bfineroot','ED_bleaf','ED_bdead','ED_balive','LEAF_MR','NPP_LEAF_SCPF', 'RH','TBOT','PBOT','QBOT','RAIN','SOILLIQ','HK','SMP','QRUNOFF','QDRAI','QRGWL','QOVER','ZWT','ZWT_PERCH','AREA_TREES','AGB','NPLANT_SCAGPFT','SEEDS_IN_LOCAL_ELEM','SEED_BANK_ELEM','SEED_GERM_ELEM','NPP_LEAF', 'M1_SCPF','M2_SCPF','M3_SCPF','M4_SCPF','M5_SCPF','M6_SCPF','M7_SCPF','M8_SCPF','MT_SCPF','CARBON_BALANCE_CANOPY_SCLS','NPLANT_CANOPY_SCPF','NPLANT_UNDERSTORY_SCPF','LAI_CANOPY_SCLS','SAI_CANOPY_SCLS','LAI_UNDERSTORY_SCLS','SAI_UNDERSTORY_SCLS','NPP_SCPF','TRIMMING_CANOPY_SCLS', 'GPP_CANOPY_SCPF','AR_CANOPY_SCPF','DDBH_CANOPY_SCPF','DDBH_UNDERSTORY_SCPF','MORTALITY_CANOPY_SCPF','NPLANT_CANOPY_SCPF','MORTALITY_UNDERSTORY_SCPF',
'FATES_TRAN_SCPF','FATES_ERRH2O_SCPF','FATES_TRAN_SCPF','FATES_SAPFLOW_SCPF','FATES_SAPFLOW_SI','FATES_ITERH1_SCPF','FATES_ITERH2_SCPF','GPP_PLANT_SCPF','G_SB_SCPF','TRIMMING_CANOPY_SCLS', 'SITE_GDD','SITE_NCHILLDAYS','SITE_NCOLDDAYS','SITE_DAYSINCE_COLDLEAFOFF','SITE_DAYSINCE_COLDLEAFON','SITE_COLD_STATUS', 'BSTOR_CANOPY_SCPF','BSTOR_UNDERSTORY_SCPF',
'FATES_ATH_SCPF', 'FATES_TTH_SCPF','FATES_STH_SCPF','FATES_LTH_SCPF','FATES_AWP_SCPF','FATES_TWP_SCPF','FATES_SWP_SCPF','FATES_LWP_SCPF', 'FATES_AFLC_SCPF','FATES_TFLC_SCPF','FATES_SFLC_SCPF','FATES_LFLC_SCPF','FATES_BTRAN_SCPF','FATES_ROOTWGT_SOILVWC_SI','FATES_ROOTWGT_SOILVWCSAT_SI', 'FATES_ROOTWGT_SOILMATPOT_SI','FATES_SOILMATPOT_SL','FATES_SOILVWC_SL','FATES_SOILVWCSAT_SL', 'FATES_ROOTUPTAKE_SI', 'FATES_ROOTUPTAKE_SL', 'FATES_ROOTUPTAKE0_SCPF','FATES_ROOTUPTAKE10_SCPF','FATES_ROOTUPTAKE50_SCPF','FATES_ROOTUPTAKE100_SCPF','H2OVEG','H2OVEG_DEAD','H2OVEG_RECRUIT', 'H2OVEG_GROWTURN_ERR','H2OVEG_PHENO_ERR','H2OVEG_HYDRO_ERR','FIRE_NESTEROV_INDEX','FIRE_ROS','FIRE_TFC_ROS','FIRE_INTENSITY','FIRE_AREA','FUEL_MOISTURE_NFSC','FIRE_FUEL_EFF_MOIST','CROWNFIREMORT_SCPF','MORTALITY_CANOPY_SCPF'
'FIRE_ROS','FIRE_TFC_ROS','FIRE_TFC_ROS_AREA_PRODUCT','FIRE_INTENSITY','FIRE_INTENSITY_AREA_PRODUCT','FIRE_AREA','FIRE_FUEL_MEF','FIRE_FUEL_BULKD', 'FIRE_FUEL_EFF_MOIST','FIRE_FUEL_SAV','SUM_FUEL','FUEL_MOISTURE_NFSC','AREA_BURNT_BY_PATCH_AGE', 'FIRE_INTENSITY_BY_PATCH_AGE','SUM_FUEL_BY_PATCH_AGE', 'BURNT_LITTER_FRAC_AREA_PRODUCT','LITTER_IN','SEED_BANK', 'SEEDS_IN'
In [1]:
# 1) import modules
import numpy as np
import pandas as pd
import datetime
from sympy import *
#
from scipy import io, integrate, linalg, signal
from scipy.sparse.linalg import eigs
from scipy.io import netcdf as nc
from scipy.optimize import curve_fit
#
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib.colors import BoundaryNorm
from matplotlib.ticker import MaxNLocator
from scipy.stats import linregress
FLUXcz1.head()
In [2]:
# read variables into array - pine
## User defined options
vars1D = ['GPP','QVEGT','EFLX_LH_TOT','RAIN','TBOT','PBOT','QBOT','RH','NPP','TLAI','AGB']
# 2D variable to explore: recruitment, seeds, DDBH,
vars2D = ['BA_SCPF','FATES_LWP_SCPF','AR_MAINT_SCPF','M2_SCPF','MT_SCPF',
'RECRUITMENT','DDBH_SCPF','NPLANT_SCPF','FATES_AWP_SCPF','NPP_SCPF',
'DDBH_CANOPY_SCPF','MORTALITY_CANOPY_SCPF','NPLANT_CANOPY_SCPF',
'DDBH_UNDERSTORY_SCPF','MORTALITY_UNDERSTORY_SCPF','NPLANT_UNDERSTORY_SCPF']
PFT=3
#define path and file surname: 1 - cz1, 2- cz2; L- co2=283, M- co2=400, H- co2=566
if PFT==1:
#for pine, use these sets
ThePFT="Pine"
pathP1L = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz1-GSWP3_1pft_Pine_Inven_T32bPn9S4_CO2Low.2021-06-29/outcome/'
fnamesurP1L = 'CZ1_PineT32bPn9S4_CO2L_'
pathP1M = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz1-GSWP3_1pft_Pine_Inven_T32bPn9S4_CO2Mid.2021-06-28/outcome/'
fnamesurP1M = 'CZ1_PineT32bPn9S4_'
pathP1H = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz1-GSWP3_1pft_Pine_Inven_T32bPn9S4_CO2High.2021-06-28/outcome/'
fnamesurP1H = 'CZ1_PineT32bPn9S4_CO2H_'
pathP2L = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz2-GSWP3_1pft_Pine_Inven_T32bPn9S4_CO2Low.2021-06-29/outcome/'
fnamesurP2L = 'CZ2_PineT32bPn9S4_CO2L_'
pathP2M = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz2-GSWP3_1pft_Pine_Inven_T32bPn9S4_CO2Mid_TP.2021-09-23/run/lnd/outcome/'
fnamesurP2M = 'CZ2_PineT32bPn9S4_CO2M_'
pathP2H = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz2-GSWP3_1pft_Pine_Inven_T32bPn9S4_CO2High.2021-06-28/run/lnd/outcome/'
fnamesurP2H = 'CZ2_PineT32bPn9S4_CO2H_'
if PFT==2:
## for oak use these sets
ThePFT="Oak"
pathP1L = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz1-GSWP3_1pft_Oak_Inven_T32bPn9S5b_CO2Low.2021-10-09/outcome/'
fnamesurP1L = 'CZ1_OakT32bPn9S5b_CO2L_'
pathP1M = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz1-GSWP3_1pft_Oak_Inven_T32bPn9S5b_CO2Mid.2021-10-08/run/lnd/outcome/'
fnamesurP1M = 'CZ1_OakT32bPn9S5b_CO2M_'
pathP1H = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz1-GSWP3_1pft_Oak_Inven_T32bPn9S5b_CO2High.2021-10-08/run/lnd/outcome/'
fnamesurP1H = 'CZ1_OakT32bPn9S5b_CO2H_'
pathP2L = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz2-GSWP3_1pft_Oak_Inven_T32bPn9S5_CO2Low.2021-10-09/outcome/'
fnamesurP2L = 'CZ2_OakT32bPn9S5b_CO2L_'
pathP2M = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz2-GSWP3_1pft_Oak_Inven_T32bPn9S5_CO2Mid.2021-10-08/run/lnd/outcome/'
fnamesurP2M = 'CZ2_OakT32bPn9S5b_CO2M_'
pathP2H = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz2-GSWP3_1pft_Oak_Inven_T32bPn9S5_CO2High.2021-10-08/run/lnd/outcome/'
fnamesurP2H = 'CZ2_OakT32bPn9S5b_CO2H_'
# pine and oak
if PFT==3:
pathP1L = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz1-GSWP3_2pft_OAK_Pine_Inven_2PFT_CO2Ls5b.2021-10-14/run/lnd/outcome/'
fnamesurP1L = 'CZ1_Pine_Oak_T32bPn9S5_'
pathP1M = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz1-GSWP3_2pft_OAK_Pine_Inven_2PFT_CO2Ms5b.2021-10-11/run/lnd/outcome/'
fnamesurP1M = 'CZ1_Pine_Oak_T32bPn9S5_'
pathP1H = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz1-GSWP3_2pft_OAK_Pine_Inven_2PFT_CO2Hs5b.2021-10-14/run/lnd/outcome/'
fnamesurP1H = 'CZ1_Pine_Oak_T32bPn9S5_'
pathP2L = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz2-GSWP3_2pft_OAK_Pine_Inven_2PFT_CO2Ls5b.2021-10-14/outcome/'
fnamesurP2L = 'CZ2_Pine_Oak_T32bPn9S5_'
pathP2M = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz2-GSWP3_2pft_OAK_Pine_Inven_2PFT_CO2Ms5b.2021-10-12/outcome/'
fnamesurP2M = 'CZ2_Pine_Oak_T32bPn9S5_'
pathP2H = '/global/scratch/junyan_ding/Models/FATES-JDgs-cz2-GSWP3_2pft_OAK_Pine_Inven_2PFT_CO2Hs5b.2021-10-14/outcome/'
fnamesurP2H = 'CZ2_Pine_Oak_T32bPn9S5_'
path=np.empty([7], dtype=object)
fnamesur=np.empty([7], dtype=object)
path[1]=pathP1L
fnamesur[1]=fnamesurP1L
path[2]=pathP1M
fnamesur[2]=fnamesurP1M
path[3]=pathP1H
fnamesur[3]=fnamesurP1H
path[4]=pathP2L
fnamesur[4]=fnamesurP2L
path[5]=pathP2M
fnamesur[5]=fnamesurP2M
path[6]=pathP2H
fnamesur[6]=fnamesurP2H
nins = 18 # total is 8 fire off + 8 fire on
baseyr=1980
styr=1980
endyr=2600
strt=(styr-baseyr)*12
endt=(endyr-baseyr)*12+12
# define 3D dimension array to store 2D variables associated with each instance
BA=np.zeros((len(path),nins+1,(endt-strt),13*3)) # basal area
GR=np.zeros((len(path),nins+1,(endt-strt),13*3)) # growth respiration
LWP=np.zeros((len(path),nins+1,(endt-strt),13*3)) # maintainence respiration
MR=np.zeros((len(path),nins+1,(endt-strt),13*3)) # maintainence respiration
MT=np.zeros((len(path),nins+1,(endt-strt),13*3)) # total mortality
M2=np.zeros((len(path),nins+1,(endt-strt),13*3)) # hydraulic mortality
DDBH=np.zeros((len(path),nins+1,(endt-strt),13*3)) # DDBH
NPLANT=np.zeros((len(path),nins+1,(endt-strt),13*3)) # NPLANT
NPPscpf=np.zeros((len(path),nins+1,(endt-strt),13*3)) # NPP by size class and pft
AWP=np.zeros((len(path),nins+1,(endt-strt),13*3)) # AWP
OTHER1scpf=np.zeros((len(path),nins+1,(endt-strt),13*3)) # other 2D variables
OTHER2scpf=np.zeros((len(path),nins+1,(endt-strt),13*3)) # other 2D variables
OTHER3scpf=np.zeros((len(path),nins+1,(endt-strt),13*3)) # other 2D variables
OTHER4scpf=np.zeros((len(path),nins+1,(endt-strt),13*3)) # other 2D variables
OTHER5scpf=np.zeros((len(path),nins+1,(endt-strt),13*3)) # other 2D variables
OTHER6scpf=np.zeros((len(path),nins+1,(endt-strt),13*3)) # other 2D variables
REC=np.zeros((len(path),nins+1,12*(endt-strt),3)) # recruitment
TLAI=np.zeros((len(path),nins+1,(endt-strt)))
AGB=np.zeros((len(path),nins+1,(endt-strt)))
GPP=np.zeros((len(path),nins+1,(endt-strt)))
NPP=np.zeros((len(path),nins+1,(endt-strt)))
QVEGT=np.zeros((len(path),nins+1,(endt-strt)))
LH=np.zeros((len(path),nins+1,(endt-strt)))
# environment variable, distinguish between cases(cid)
RAIN=np.zeros((len(path),(endt-strt))) # rain
TBOT=np.zeros((len(path),(endt-strt))) # temperature
RYear=np.zeros((len(path),(endt-strt)))
RMonth=np.zeros((len(path),(endt-strt)))
MaxT=0 #maximum length of time
#GR = np.array[[[]]]
#AR = np.array[[[]]]
cases=np.array([1,2,3,4,5,6])
# cid refer to the treatment: 2 cz sites X 3 CO2 level
for cid in cases:
for ins in range (1,nins+1):
Theins = format(ins, '04')
filenm = path[cid] + fnamesur[cid] + Theins +'.nc'
print(cid, path[cid], fnamesur[cid])
print(filenm)
fin = nc.netcdf_file(filenm)
### read general variables if this is the first instance
Thetime = np.array(fin.variables['time'][:])
for iTime in range(0,int(len(Thetime)/12)):
for iMon in range (1,13):
RYear[cid][iTime*12+iMon-1] = styr + iTime # convert fates time fueld to year
RMonth[cid][(iTime*12+iMon-1)] = iMon # conver fates time field to month
# read size class, pft number, and environment variables when it is the first instance
if ins == 1 :
Thesizebin = np.array(fin.variables['fates_levscls'][:])
Npft = len(np.squeeze(np.array(fin.variables['fates_levpft'][:])))
print(Npft)
BAtmp = np.squeeze(np.array(fin.variables[vars2D[0]][0:len(Thetime)]))
RECtmp = np.squeeze(np.array(fin.variables[vars2D[5]][0:len(Thetime)]))
MTtmp = np.squeeze(np.array(fin.variables[vars2D[4]][0:len(Thetime)]))
M2tmp = np.squeeze(np.array(fin.variables[vars2D[3]][0:len(Thetime)]))
DDBHtmp = np.squeeze(np.array(fin.variables[vars2D[6]][0:len(Thetime)]))
NPLANTtmp = np.squeeze(np.array(fin.variables[vars2D[7]][0:len(Thetime)]))
AWPtmp = np.squeeze(np.array(fin.variables[vars2D[8]][0:len(Thetime)]))
LWPtmp = np.squeeze(np.array(fin.variables[vars2D[1]][0:len(Thetime)]))
NPPscpftmp = np.squeeze(np.array(fin.variables[vars2D[9]][0:len(Thetime)]))
OTHER1scpftmp = np.squeeze(np.array(fin.variables[vars2D[10]][0:len(Thetime)]))
OTHER2scpftmp = np.squeeze(np.array(fin.variables[vars2D[11]][0:len(Thetime)]))
OTHER3scpftmp = np.squeeze(np.array(fin.variables[vars2D[12]][0:len(Thetime)]))
OTHER4scpftmp = np.squeeze(np.array(fin.variables[vars2D[13]][0:len(Thetime)]))
OTHER5scpftmp = np.squeeze(np.array(fin.variables[vars2D[14]][0:len(Thetime)]))
OTHER6scpftmp = np.squeeze(np.array(fin.variables[vars2D[15]][0:len(Thetime)]))
BA[cid][ins][0:len(Thetime)][:] = BAtmp[0:len(Thetime)][:]
REC[cid][ins][0:len(Thetime)][:] = RECtmp[0:len(Thetime)][:]
MT[cid][ins][0:len(Thetime)][:] = MTtmp[0:len(Thetime)][:]
DDBH[cid][ins][0:len(Thetime)][:] = DDBHtmp[0:len(Thetime)][:]
NPLANT[cid][ins][0:len(Thetime)][:] = NPLANTtmp[0:len(Thetime)][:]
AWP[cid][ins][0:len(Thetime)][:] = AWPtmp[0:len(Thetime)][:]
LWP[cid][ins][0:len(Thetime)][:] = LWPtmp[0:len(Thetime)][:]
NPPscpf[cid][ins][0:len(Thetime)][:] = NPPscpftmp[0:len(Thetime)][:]
OTHER1scpf[cid][ins][0:len(Thetime)][:] = OTHER1scpftmp[0:len(Thetime)][:]
OTHER2scpf[cid][ins][0:len(Thetime)][:] = OTHER2scpftmp[0:len(Thetime)][:]
OTHER3scpf[cid][ins][0:len(Thetime)][:] = OTHER3scpftmp[0:len(Thetime)][:]
OTHER4scpf[cid][ins][0:len(Thetime)][:] = OTHER4scpftmp[0:len(Thetime)][:]
OTHER5scpf[cid][ins][0:len(Thetime)][:] = OTHER5scpftmp[0:len(Thetime)][:]
OTHER6scpf[cid][ins][0:len(Thetime)][:] = OTHER6scpftmp[0:len(Thetime)][:]
TLAI[cid][ins][0:len(Thetime)] = np.squeeze(np.array(fin.variables[vars1D[9]][0:len(Thetime)]))
AGB[cid][ins][0:len(Thetime)] = np.squeeze(np.array(fin.variables[vars1D[10]][0:len(Thetime)]))
GPP[cid][ins][0:len(Thetime)] = np.squeeze(np.array(fin.variables[vars1D[0]][0:len(Thetime)]))
NPP[cid][ins][0:len(Thetime)] = np.squeeze(np.array(fin.variables[vars1D[8]][0:len(Thetime)]))
QVEGT[cid][ins][0:len(Thetime)] = np.squeeze(np.array(fin.variables[vars1D[1]][0:len(Thetime)]))
LH[cid][ins][0:len(Thetime)] = np.squeeze(np.array(fin.variables[vars1D[2]][0:len(Thetime)]))
RAIN[cid][0:len(Thetime)] = np.squeeze(np.array(fin.variables[vars1D[3]][0:len(Thetime)]))
TBOT[cid][0:len(Thetime)] = np.squeeze(np.array(fin.variables[vars1D[4]][0:len(Thetime)]))
fin.close()
BApftsum=np.zeros((len(path),nins+1,(endt-strt),3))
BAtotal=np.zeros((len(path),nins+1,(endt-strt)))
BApctPine=np.zeros((len(path),nins+1,(endt-strt)))
MTbapftsum=np.zeros((len(path),nins+1,(endt-strt),3))
DBAtotal=np.zeros((len(path),nins+1,(endt-strt)))
MTbatotal=np.zeros((len(path),nins+1,(endt-strt)))
M2batotal=np.zeros((len(path),nins+1,(endt-strt)))
# TOTAL growth respiration by pft
GRpftsum=np.zeros((len(path),nins+1,(endt-strt),3))
GRtotal=np.zeros((len(path),nins+1,(endt-strt)))
GRpctPine=np.zeros((len(path),nins+1,(endt-strt)))
# total maintenance respiration by pft
MRpftsum=np.zeros((len(path),nins+1,(endt-strt),3))
MRtotal=np.zeros((len(path),nins+1,(endt-strt)))
MRpctPine=np.zeros((len(path),nins+1,(endt-strt)))
#ADJUST AWP and MTba for nan by NPLNAT
AWPadj=(AWP/NPLANT)*NPLANT
LWPadj=(LWP/NPLANT)*NPLANT
DDBHadj=DDBH/NPLANT
MTadj=(MT/NPLANT)*NPLANT
NPPscpfadj=NPPscpf/NPLANT
In [3]:
# adjust for nplant for other 2D variables,
# others:DDBH canopy, Mortality canopy, nplant canopy
# others:DDBH understory, mortality understory, nplant understory
OTHER1scpfadj=OTHER1scpf/OTHER3scpf # DDBH per tree canopy
OTHER2scpfadj=OTHER2scpf/OTHER3scpf # mortality rate canopy
OTHER4scpfadj=OTHER4scpf/OTHER6scpf # DDBY per tree understory
OTHER5scpfadj=OTHER5scpf/OTHER6scpf # mortality rate understory
NPLANT_CAN=OTHER3scpf
len(Thetime)/12¶
print(NPLANTtmp.shape)¶
ArrayTst=np.zeros((1800,39))¶
print("ArrayTst")¶
print(ArrayTst[:,1])¶
print(ArrayTst[1][:])¶
MTadj=(MT/NPLANT)*NPLANT print (Thetime[1200], 5400/12)
print(RYear[3][0:100])¶
print(GPP[3][1][0:100] )¶
print(NPLANTtmp.shape) print("NPLANTtmp[3,:]:",NPLANTtmp[3,:]) print("MT[5][23][400:500,5]",MT[5][23][400:500,5]) print("MTtmp[:,5]:",MTtmp[:,5]) print("AWP") print(AWP[5][1][1:15][0:12]) print("AWPadj")
AWPadj=(AWP/NPLANT)*NPLANT¶
print(AWPadj[5][1][3][0:12])
NPLANT[5][21][0:300][8]¶
AWPtest=np.zeros((len(path),nins+1,(endt-strt),13*3)) AWPtest[5][1][3][0:12]=AWPadj[5][1][3][0:12] print("====") print(np.nanmean(AWPtest[5][1][3][0:12])) print(np.mean(AWPtest[5][1][3][0:12]))
In [4]:
DBH=np.zeros(len(Thesizebin)*2)
CIRC=np.zeros(len(Thesizebin)*2)
for i in range(0,len(Thesizebin)-1):
DBH[i]=(Thesizebin[i]+Thesizebin[i+1])/2
DBH[i+13]=(Thesizebin[i]+Thesizebin[i+1])/2
# diameter
DBH[12]=Thesizebin[12]+5
DBH[12+13]=Thesizebin[12]+5
BAchr=(3.14*(DBH/2)**2)/10000 # 2*3.14*R
#GPP[15]
#BAchr=BA/NPLANT # can't use this method for canopy and understory
#DBHchr=2*(BAchr/3.14)**0.5 # can't use this method for canopy and understory
BAtree=(3.14*(DBH/2)**2)/10000 # basal area of the individual
DDBH_CANscpfadj=OTHER1scpfadj
MT_CANscpfadj=OTHER2scpfadj
DDBH_UNDscpfadj=OTHER4scpfadj
MT_UNDscpfadj=OTHER5scpfadj
# DBH is in cm, BA is in m^2, so divide by 10000 to convert cm^2 to m^2
def CAL_DBA(DBH, DDBH):
DBA=(3.14*(DBH/2+DDBH/2)**2 - 3.14*(DBH/2)**2)/10000
return DBA
In [31]:
print(DBH[:])
print((3.14*(DBH[:]/2)**2)/10000)
#print(Thesizebin)
print("BAchr[5][23][1:150,6]",BAchr[5][23][1:150,6])
print("DBHchr[5][23][1:150,5]",DBHchr[5][23][1:150,5])
print("underestimation of DBA (cm^2) per tree for a 39 cm DBH with DDBH=0.5cm:")
print((3.14*(39/2+0.25)**2-3.14*(39/2)**2)-(3.14*(35.2/2+0.25)**2-3.14*(35.2/2)**2))
print("actuall DBA ", 3.14*(39/2+0.25)**2-3.14*(39/2)**2)
print("underestimation of DBA (cm^2) per tree for a 39 cm DBH with DDBH=0.8cm:")
print((3.14*(39/2+0.4)**2-3.14*(39/2)**2)-(3.14*(35.2/2+0.4)**2-3.14*(35.2/2)**2))
print("actuall DBA (m^2)", (3.14*(39/2+0.4)**2-3.14*(39/2)**2)/10000)
In [19]:
# 1) calculate the total BA and MR, AR, GPP, MT, AWP of each pft
# 2) calculate 33 years statistic:mean, variation, SD ,of BA, MTba : 100 to 133 year by case and instances
# 3) calculate annual mean BA, MTba, AWP
# 4) calculate DDBH Recruitment to MTba ratio
# define time range for statistic
# DIMA is the diameter of the cohort, DBA=DDBH*DIMA of each individual in that size class
# total DBA of a cohort is DBAtotal=TDBA*NPLANTS
endyr=150
styr=endyr-33
# zero BA, MTba statis variables
DBA=np.zeros((len(path),nins+1,(endt-strt),13*3)) # DBA by month by cohort
DBA_CAN=np.zeros((len(path),nins+1,(endt-strt),13*3)) # DBA by month by cohort
DBA_UND=np.zeros((len(path),nins+1,(endt-strt),13*3)) # DBA by month by cohort
MTba_CAN=np.zeros((len(path),nins+1,(endt-strt),13*3)) # DBA by month by cohort
MTba_UND=np.zeros((len(path),nins+1,(endt-strt),13*3)) # DBA by month by cohort
DBAtotal=np.zeros((len(path),nins+1,(endt-strt)))
BAmean=np.zeros((7,nins+1))
BAstd=np.zeros((7,nins+1))
MTbamean=np.zeros((7,nins+1))
MTbastd=np.zeros((7,nins+1))
M2bamean=np.zeros((7,nins+1))
M2bastd=np.zeros((7,nins+1))
DBAmean=np.zeros((7,nins+1)) #increase in BA area as converted from DDBH m^2
DBAstd=np.zeros((7,nins+1)) #
# annual mean
DDBHadjyrmean=np.zeros((7,nins+1,endyr+2,13*3))
DBAyrmean=np.zeros((7,nins+1,endyr+2,13*3))
BAyrmean=np.zeros((7,nins+1,endyr+2,13*3))
MTyrmean=np.zeros((7,nins+1,endyr+2,13*3))
MTbayrmean=np.zeros((7,nins+1,endyr+2,13*3))
NPLANTyrmean=np.zeros((7,nins+1,endyr+2,13*3))
NPLANT_CANyrmean=np.zeros((7,nins+1,endyr+2,13*3))
DDBH_CANadjyrmean=np.zeros((7,nins+1,endyr+2,13*3))
DBA_CANadjyrmean=np.zeros((7,nins+1,endyr+2,13*3))
MTba_CANadjyrmean=np.zeros((7,nins+1,endyr+2,13*3))
DDBH_UNDadjyrmean=np.zeros((7,nins+1,endyr+2,13*3))
DBA_UNDadjyrmean=np.zeros((7,nins+1,endyr+2,13*3))
MTba_UNDadjyrmean=np.zeros((7,nins+1,endyr+2,13*3))
AWPyrmean=np.zeros((7,nins+1,endyr+2,13*3))
LWPyrmean=np.zeros((7,nins+1,endyr+2,13*3))
TLAIyrmean=np.zeros((7,nins+1,endyr+2))
AGByrmean=np.zeros((7,nins+1,endyr+2))
GPPyrmean=np.zeros((7,nins+1,endyr+2))
QVEGTyrmean=np.zeros((7,nins+1,endyr+2))
TBAyrmean=np.zeros((7,nins+1,endyr+2))
RAINyrmean=np.zeros((7,endyr+2)) # annual mean precipitation rate mm/year
TBOTyrmean=np.zeros((7,endyr+2)) # annual mean temperature
# BA weighted DDBHmean & std, MTmean & std,
DBAmeanByBA=np.zeros((7,nins+1))
DBAstdByBA=np.zeros((7,nins+1))
MTbameanByBA=np.zeros((7,nins+1))
MTbastdByBA=np.zeros((7,nins+1))
# by pft, population level
LWPpftyrmean=np.zeros((len(path),nins+1,endyr+2,4))
AWPpftyrmean=np.zeros((len(path),nins+1,endyr+2,4))
TBApftyrmean=np.zeros((len(path),nins+1,endyr+2,4))
RECbamean=np.zeros((7,nins+1))
BArec=3.14*((2.5/100)/2)**2 # basal of the individual of the recruitment, dbh=2.5cm
SperYr=365*24*60*60
# convert MT to basal area
MTba=BA*MT
M2ba=BA*M2
# Individual level
DDBHadj=DDBH/NPLANT # DDBH in fates history file is calculated as the sum of all the individuals in that size class
MTbaadj=MTba/NPLANT
MRch=MR*(BA/NPLANT)
# calculating
for cid in cases:
for theins in range (1,nins+1):
#print ('theins:', theins)
for Ctime in range (1,(endt-1)):
for Thepft in range (1,3):
BApftsum[cid][theins][Ctime][Thepft]=np.sum(BA[cid][theins][Ctime,((Thepft-1)*13):(Thepft-1)*13+12])
MRpftsum[cid][theins][Ctime][Thepft]=np.sum(MR[cid][theins][Ctime,((Thepft-1)*13):(Thepft-1)*13+12])
GRpftsum[cid][theins][Ctime][Thepft]=np.sum(GR[cid][theins][Ctime,((Thepft-1)*13):(Thepft-1)*13+12])
MTbapftsum[cid][theins][Ctime][Thepft]=np.sum(MTba[cid][theins][Ctime,((Thepft-1)*13):(Thepft-1)*13+12])
#print ('thepft:',Thepft)
# calculate total BA, DBA, MT for each time step of the frid
# since this is one PFT runs, the grid total also is the total of the PFT
BAtotal[cid][theins][Ctime]=np.sum((BA[cid][theins][Ctime,0:25]))
# unit of DBAtotal is m^2, converted from DDBH(cm) and DIMA (cm),
# (DIMA+0.5DDBH)*NPLANT
DBA[cid][theins][Ctime][0:25]=(3.14*(DDBHadj[cid][theins][Ctime][0:25]/2+DBH[0:25]/2)**2-3.14*(DBH[0:25]/2)**2)/10000
DBA_CAN[cid][theins][Ctime][0:25]=CAL_DBA(BAtree[0:25],DDBH_CANscpfadj[cid][theins][Ctime][0:25])
DBA_UND[cid][theins][Ctime][0:25]=CAL_DBA(BAtree[0:25],DDBH_UNDscpfadj[cid][theins][Ctime][0:25])
MTba_CAN[cid][theins][Ctime][0:25]=MT_CANscpfadj[cid][theins][Ctime][0:25]*BAtree[0:25]
MTba_UND[cid][theins][Ctime][0:25]=MT_UNDscpfadj[cid][theins][Ctime][0:25]*BAtree[0:25]
DBAtotal[cid][theins][Ctime]=np.nansum((3.14*(DDBHadj[cid][theins][Ctime][0:25]/2+DBH[0:25]/2)**2-3.14*(DBH[0:25]/2)**2))/10000
BApctPine[cid][theins][Ctime]=BApftsum[cid][theins][Ctime][1]/BAtotal[cid][theins][Ctime]
MRtotal[cid][theins][Ctime]=np.sum((MR[cid][theins][Ctime][0:25]))
MRpctPine[cid][theins][Ctime]=MRpftsum[cid][theins][Ctime][1]/MRtotal[cid][theins][Ctime]
GRpctPine[cid][theins][Ctime]=GRpftsum[cid][theins][Ctime][1]/GRtotal[cid][theins][Ctime]
MTbatotal[cid][theins][Ctime]=np.sum((MTba[cid][theins][Ctime][0:25]))
M2batotal[cid][theins][Ctime]=np.sum((M2ba[cid][theins][Ctime][0:25]))
# CALCULATE annual average of DBA, MT, AWPadj of each cohort, ingore the last year
for Ytime in range (1, (endyr-1)):
# I think there is a bug in the python as to check the array length. If I put
#
#startTime=(Ytime-1)*12+1
#endTime=(Ytime-1)*12+15
# set the value to be 10 12+10
strMon=10
endMon=11+10
#
for chid in range (0,26):
#print("year", Ytime,"startime:", startTime, "endTime:", endTime)
#print(chid,Ctime,DDBH[cid][theins][startTime:endTime][chid])
DDBHadjyrmean[cid][theins][Ytime][chid]=np.nanmean(DDBHadj[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
BAyrmean[cid][theins][Ytime][chid]= np.nanmean(BA[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
DBAyrmean[cid][theins][Ytime][chid]=np.nanmean(DBA[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
MTbayrmean[cid][theins][Ytime][chid]=np.nanmean(MTbaadj[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
NPLANTyrmean[cid][theins][Ytime][chid]=np.nanmean(NPLANT[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
DBA_CANadjyrmean[cid][theins][Ytime][chid]=np.nanmean(DBA_CAN[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
MTba_CANadjyrmean[cid][theins][Ytime][chid]=np.nanmean(MTba_CAN[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
DBA_UNDadjyrmean[cid][theins][Ytime][chid]=np.nanmean(DBA_UND[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
MTba_UNDadjyrmean[cid][theins][Ytime][chid]=np.nanmean(MTba_UND[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
DDBH_CANadjyrmean[cid][theins][Ytime][chid]=np.nanmean(DDBH_CANscpfadj[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
DDBH_UNDadjyrmean[cid][theins][Ytime][chid]=np.nanmean(DDBH_UNDscpfadj[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
NPLANT_CANyrmean[cid][theins][Ytime][chid]=np.nanmean(NPLANT_CAN[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
LWPyrmean[cid][theins][Ytime][chid]=np.nanmean(LWPadj[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
AWPyrmean[cid][theins][Ytime][chid]=np.nanmean(AWPadj[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
MTyrmean[cid][theins][Ytime][chid]=np.nanmean(MTadj[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
# calculater for pft1 obnly
for Thepft in range (1,3):
LWPpftyrmean[cid][theins][Ytime][Thepft]=np.nanmean(LWPadj[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),((Thepft-1)*13):(Thepft-1)*13+12])
AWPpftyrmean[cid][theins][Ytime][Thepft]=np.nanmean(AWPadj[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),((Thepft-1)*13):(Thepft-1)*13+12])
TBApftyrmean[cid][theins][Ytime][Thepft]=np.nanmean(BApftsum[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),Thepft])
TBAyrmean[cid][theins][Ytime]=np.nanmean(BAtotal[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon)])
RAINyrmean[cid][Ytime]=np.mean(RAIN[cid][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon)])*SperYr
TBOTyrmean[cid][Ytime]=np.mean(TBOT[cid][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon)])
TLAIyrmean[cid][theins][Ytime]=np.nanmean(TLAI[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon)])
AGByrmean[cid][theins][Ytime]=np.nanmean(AGB[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon)])
GPPyrmean[cid][theins][Ytime]=SperYr*np.nanmean(GPP[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon)])
QVEGTyrmean[cid][theins][Ytime]=SperYr*np.nanmean(QVEGT[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon)])
# calculate mean and std for the last 33 years
#print(REC[cid][theins][styr*12:endyr*12][0])
RECbamean[cid][theins]=np.mean(REC[cid][theins][styr*12:endyr*12,0])*BArec
BAmean[cid][theins]=np.mean(BAtotal[cid][theins][styr*12:endyr*12])
BAstd[cid][theins]=np.std(BAtotal[cid][theins][styr*12:endyr*12])
MTbamean[cid][theins]=np.mean(MTbatotal[cid][theins][styr*12:endyr*12])
MTbastd[cid][theins]=np.std(MTbatotal[cid][theins][styr*12:endyr*12])
DBAmean[cid][theins]=np.mean(DBAtotal[cid][theins][styr*12:endyr*12])
DBAstd[cid][theins]=np.std(DBAtotal[cid][theins][styr*12:endyr*12])
#BA weighted
MTbameanByBA[cid][theins]=np.mean((MTbatotal[cid][theins][styr*12:endyr*12]/BAtotal[cid][theins][styr*12:endyr*12]))
MTbastdByBA[cid][theins]=np.std((MTbatotal[cid][theins][styr*12:endyr*12]/BAtotal[cid][theins][styr*12:endyr*12]))
DBAmeanByBA[cid][theins]=np.mean((DBAtotal[cid][theins][styr*12:endyr*12]/BAtotal[cid][theins][styr*12:endyr*12]))
DBAstdByBA[cid][theins]=np.std((DBAtotal[cid][theins][styr*12:endyr*12]/BAtotal[cid][theins][styr*12:endyr*12]))
# compute the difference in BAmean,std; DDBHmean,DDBHstd; MTmean, std between low-medium,
# high - medium co2 by per ppm change of c02
DifCO2_low_med=386-283
DifCO2_high_med=566-386
In [ ]:
# plot lwp vs awp
In [27]:
# compute DDBH and MTba of size class of canopy tree
# DBA_CANadjyrmean[cid][theins][Ytime][chid]=np.nanmean(DBA_CAN[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
# MTba_CANadjyrmean[cid][theins][Ytime][chid]=np.nanmean(MTba_CAN[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
# DBA_UNDadjyrmean[cid][theins][Ytime][chid]=np.nanmean(DBA_UND[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
# MTba_UNDadjyrmean[cid][theins][Ytime][chid]=np.nanmean(MTba_UND[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
# DDBH_CANadjyrmean[cid][theins][Ytime][chid]=np.nanmean(DDBH_CANscpfadj[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
# DDBH_UNDadjyrmean[cid][theins][Ytime][chid]=np.nanmean(DDBH_UNDscpfadj[cid][theins][((Ytime-1)*12+strMon):((Ytime-1)*12+endMon),chid])
sets=np.array([1,3,6])
Thechid=0 # the cohort size id
# create an empty datafram
# define column names
recid=0 # table key index
column_names = ["cid", "instid","PFT", "chid","DBAmn","MTbamn"]
ChStc = pd.DataFrame(columns = column_names)
for cid in cases:
for theins in sets:
for chid in range (0,26):
#print(recid)
cname="DBA_c"+str(cid)+"s"+str(theins)
DBAmn=10000*np.nanmean(DBA_CANadjyrmean[cid][theins][70:149,chid])
MTbamn=10000*np.nanmean(MTba_CANadjyrmean[cid][theins][70:149,chid])
if chid<=12:
thepft=1
else:
thepft=2
#print(MTbamn)
ChStc.loc[recid]=[cid,theins,thepft,chid,DBAmn,MTbamn]
recid=recid+1
recid=0 # table key index
column_names = ["cid", "instid","PFT", "chid","DBAmn","MTbamn"]
ChStcUND = pd.DataFrame(columns = column_names)
for cid in cases:
for theins in sets:
for chid in range (0,26):
#print(recid)
cname="DBA_c"+str(cid)+"s"+str(theins)
DBAmn=10000*np.nanmean(DBA_UNDadjyrmean[cid][theins][70:149,chid])
MTbamn=10000*np.nanmean(MTba_UNDadjyrmean[cid][theins][70:149,chid])
if chid<=12:
thepft=1
else:
thepft=2
#print(MTbamn)
ChStcUND.loc[recid]=[cid,theins,thepft,chid,DBAmn,MTbamn]
recid=recid+1
In [24]:
print(ChStc)
In [25]:
# plot mean DBA vs MTba of the canopy level
for theins in [3]:
ChStcSel=ChStc[(ChStc.cid.isin([1,2,3,4,5,6])) & (ChStc.instid.isin([theins]))]
ChStcSel0=ChStcSel[ChStcSel.chid.isin([0,13])]
ChStcSel1=ChStcSel[ChStcSel.chid.isin([1,14])]
ChStcSel5=ChStcSel[ChStcSel.chid.isin([5,18])]
ChStcSel8=ChStcSel[ChStcSel.chid.isin([8,21])]
ChStcSel11=ChStcSel[ChStcSel.chid.isin([11,24])]
ax1 = ChStcSel0.plot.scatter(x='MTbamn',y='DBAmn',s=(ChStcSel0["cid"]+2)**2, c='PFT',colormap='Spectral')
ax1.set_title("size 2.5cm DBH ins:"+str(theins))
ax1.set_xlabel("MTba")
ax2 = ChStcSel1.plot.scatter(x='MTbamn',y='DBAmn',s=(ChStcSel1["cid"]+2)**2, c='PFT',colormap='Spectral')
ax2.set_title("size 7.5cm DBH ins:"+str(theins))
ax3 = ChStcSel1.plot.scatter(x='MTbamn',y='DBAmn',s=(ChStcSel5["cid"]+2)**2, c='PFT',colormap='Spectral')
ax3.set_title("size 35cm DBH ins:"+str(theins))
ax4 = ChStcSel8.plot.scatter(x='MTbamn',y='DBAmn',s=(ChStcSel8["cid"]+2)**2, c='PFT',colormap='Spectral')
ax4.set_title("size 65cm DBH ins:"+str(theins))
ax5 = ChStcSel11.plot.scatter(x='MTbamn',y='DBAmn',s=(ChStcSel11["cid"]+2)**2, c='PFT',colormap='Spectral')
ax5.set_title("size 85cm DBH ins:"+str(theins))
In [28]:
# plot mean DBA vs MTba of understory together
for theins in [3]:
ChStcSel=ChStcUND[(ChStc.cid.isin([1,2,3,4,5,6])) & (ChStc.instid.isin([theins]))]
ChStcSel0=ChStcSel[ChStcSel.chid.isin([0,13])]
ChStcSel1=ChStcSel[ChStcSel.chid.isin([1,14])]
ChStcSel5=ChStcSel[ChStcSel.chid.isin([5,18])]
ChStcSel8=ChStcSel[ChStcSel.chid.isin([8,21])]
ChStcSel11=ChStcSel[ChStcSel.chid.isin([11,24])]
ax1 = ChStcSel0.plot.scatter(x='MTbamn',y='DBAmn',s=(ChStcSel0["cid"]+2)**2, c='PFT',colormap='Spectral')
ax1.set_title("size 2.5cm DBH ins:"+str(theins))
ax1.set_xlabel("MTba")
ax2 = ChStcSel1.plot.scatter(x='MTbamn',y='DBAmn',s=(ChStcSel1["cid"]+2)**2, c='PFT',colormap='Spectral')
ax2.set_title("size 7.5cm DBH ins:"+str(theins))
ax3 = ChStcSel1.plot.scatter(x='MTbamn',y='DBAmn',s=(ChStcSel5["cid"]+2)**2, c='PFT',colormap='Spectral')
ax3.set_title("size 35cm DBH ins:"+str(theins))
ax4 = ChStcSel8.plot.scatter(x='MTbamn',y='DBAmn',s=(ChStcSel8["cid"]+2)**2, c='PFT',colormap='Spectral')
ax4.set_title("size 65cm DBH ins:"+str(theins))
ax5 = ChStcSel11.plot.scatter(x='MTbamn',y='DBAmn',s=(ChStcSel11["cid"]+2)**2, c='PFT',colormap='Spectral')
ax5.set_title("size 85cm DBH ins:"+str(theins))
In [19]:
# plot mean DBA vs MTba
for theins in [1,3,6]:
ChStcSel=ChStc[(ChStc.cid.isin([1,2,3,4,5,6])) & (ChStc.instid.isin([theins]))]
ChStcSel0=ChStcSel[ChStcSel['chid']==0 ]
ChStcSel1=ChStcSel[ChStcSel['chid']==1 ]
ChStcSel5=ChStcSel[ChStcSel['chid']==2 ]
ChStcSel8=ChStcSel[ChStcSel['chid']==8 ]
ChStcSel13=ChStcSel[ChStcSel['chid']==13 ]
ChStcSel14=ChStcSel[ChStcSel['chid']==14 ]
ChStcSel19=ChStcSel[ChStcSel['chid']==15 ]
ChStcSel21=ChStcSel[ChStcSel['chid']==21 ]
ax1 = ChStcSel0.plot.scatter(x='MTbamn',y='DBAmn',s=ChStcSel0["cid"]*10, c='cid',colormap='viridis')
ax1.set_title("size 2.5cm DBH pine")
ax1.set_xlabel("MTba")
ax2 = ChStcSel1.plot.scatter(x='MTbamn',y='DBAmn',s=ChStcSel1["cid"]*10, c='cid',colormap='viridis')
ax2.set_title("size 7.5cm DBH pine")
ax3 = ChStcSel8.plot.scatter(x='MTbamn',y='DBAmn',s=ChStcSel8["cid"]*10, c='cid',colormap='viridis')
ax3.set_title("size 65cm DBH pine")
ax4 = ChStcSel13.plot.scatter(x='MTbamn',y='DBAmn',s=ChStcSel13["cid"]*10, c='cid',colormap='viridis')
ax4.set_title("size 2.5cm DBH oak")
ax4.set_xlabel("MTba")
ax5 = ChStcSel14.plot.scatter(x='MTbamn',y='DBAmn',s=ChStcSel14["cid"]*10, c='cid',colormap='viridis')
ax5.set_title("size 7.5cm DBH oak")
ax6 = ChStcSel21.plot.scatter(x='MTbamn',y='DBAmn',s=ChStcSel21["cid"]*10, c='cid',colormap='viridis')
ax6.set_title("size 65cm DBH oak")
ax4.legend()
In [9]:
#plot histogram or boxplot of TBA, annual GPP and ET for each instance
# TBAyrmean[cid][theins][Ytime]
Bchid=12
for cid in cases:
fig, bp = plt.subplots()
bpdata0=TBAyrmean[cid][sets[0]][1:148]
bpdata1=TBAyrmean[cid][sets[1]][1:148]
bpdata2=TBAyrmean[cid][sets[2]][1:148]
bpdata3=TBAyrmean[cid][sets[3]][1:148]
bp.boxplot([bpdata0,bpdata1,bpdata2,bpdata3])
bp.set_ylabel("TBA (m^2/ha)")
plt.title("Total basal area (case"+str(cid)+")")
# plot size distribution of canopy level as a color map
for cid in cases:
fig, bps = plt.subplots()
bpdata0=NPLANT_CANyrmean[cid][sets[0]][100:148,Bchid+pftadj]
bpdata1=NPLANT_CANyrmean[cid][sets[1]][100:148,Bchid+pftadj]
bpdata2=NPLANT_CANyrmean[cid][sets[2]][100:148,Bchid+pftadj]
bpdata3=NPLANT_CANyrmean[cid][sets[3]][100:148,Bchid+pftadj]
bps.boxplot([bpdata0,bpdata1,bpdata2,bpdata3])
plt.title("number large trees (DBH 85 cm),(case"+str(cid)+")")
bps.set_ylabel("N trees (n/ha)")
In [11]:
# GPP, QVEGT vs. TBA
Thestyr=1
Theendyr=149
for cid in cases:
fig,GPPpt=plt.subplots(nrows=2,ncols=2,figsize=(11,11))
insid=0
for i in range(0,2):
for j in range(0,2):
# plot DBA vs MTba
GPPpt[i][j].scatter(TBAyrmean[cid][sets[insid]][Thestyr:Theendyr],GPPyrmean[cid][sets[insid]][Thestyr:Theendyr]) # ,marker=MK[1:54]
# BAvsMT[i][j].scatter(AWPyrmean[Thecid][inset[insid]][Thestyr:Theendyr,Thechid-max(i,j)], MTbayrmean[cid][inset[insid]][Thestyr:Theendyr,Thechid-max(i,j)]) # ,marker=MK[1:54]
GPPpt[i][j].set_xlabel('TBA (cm^2/ha)', labelpad = 5)
GPPpt[i][j].set_ylabel('GPP (gC/s/ha)', labelpad = 5)
GPPpt[i][j].set_title("GPP vs TBA: c"+str(cid)+"ins"+str(sets[insid]))
#BAvsMT[i][j].set_ylim([0, 0.00002])
#BAvsMT[i][j].set_xlim([0, 0.00002])
insid=insid+1
fig,QVEGTpt=plt.subplots(nrows=2,ncols=2,figsize=(11,11))
insid=0
for i in range(0,2):
for j in range(0,2):
# plot DBA vs MTba
QVEGTpt[i][j].scatter(TBAyrmean[cid][sets[insid]][Thestyr:Theendyr],QVEGTyrmean[cid][sets[insid]][Thestyr:Theendyr]) # ,marker=MK[1:54]
# BAvsMT[i][j].scatter(AWPyrmean[Thecid][inset[insid]][Thestyr:Theendyr,Thechid-max(i,j)], MTbayrmean[cid][inset[insid]][Thestyr:Theendyr,Thechid-max(i,j)]) # ,marker=MK[1:54]
QVEGTpt[i][j].set_xlabel('TBA (m^2/ha)', labelpad = 5)
QVEGTpt[i][j].set_ylabel('QVEGT (m/s/ha)', labelpad = 5)
QVEGTpt[i][j].set_title("transpiration vs TBA: c"+str(cid)+"ins"+str(sets[insid]))
#BAvsMT[i][j].set_ylim([0, 0.00002])
#BAvsMT[i][j].set_xlim([0, 0.00002])
insid=insid+1
In [ ]:
## plot total basal area of two pft by month
## dinguish pft by line type, instance by color
cases=[1,2,3,4,5,6]
sets=np.array([1,3,6])
## assign line type, and line color
lc = np.array(['blue','red','green','purple','orange','black','olive','brown','peachpuff']) # line color
ls = np.array(['-','--',])
## plot
## setup plots to make
lid1=0
lid2=0
fig,(pBA, Bt)=plt.subplots(nrows=2,ncols=1,figsize=(20,15))
for cid in cases:
if inst <=Oneins:
ils = 0
ilc = lid1
lid1=lid1+1
else:
ils = 1
ilc = lid2
lid2=lid2+1
pBA.plot(Thetime[1:len(Thetime)], BApctPine[inst][1:(len(Thetime))], linestyle=ls[ils], color=lc[ilc], linewidth=1.5)
#pBA.plot(Thetime[1:len(Thetime)],[0.3,0.3],'b-')
pBA.set_xlabel("day since 1980-01-01")
pBA.set_ylabel("BApine/BAtotal")
pBA.legend(cases, title="Instance", loc='upper left')
lid1=0
lid2=0
for inst in cases:
if inst <=Oneins:
ils = 0
ilc = lid1
lid1=lid1+1
else:
ils = 1
ilc = lid2
lid2=lid2+1
Bt.plot(Thetime[1:len(Thetime)], BAtotal[inst][1:(len(Thetime))], linestyle=ls[ils], color=lc[ilc], linewidth=1.5)
Bt.set_xlabel("day since 1980-01-01")
Bt.set_ylabel("BAtotal (m^2/ha)")
Bt.legend(cases, title="Instance",loc='upper left')
#plt.show
#print("notes: large fluctuation are because of the storagy cushion is too low - 2, so during a cold year, most of oak are died from carbon starvation, while during normal year, they have more carbon to grow")
In [13]:
# plot 2D var as map by cohort x time
pcid=5
pinst=35 #5 11 29 35
pft=PFT
pftadj=pft*13-13
Z=np.transpose(DDBHadjyrmean[pcid][pinst][10:150,0+pftadj:12+pftadj])
fig, figplt = plt.subplots(nrows=1,ncols=1, figsize=(20,40))
figplt.imshow(Z)
figplt.set_title("DDBH per year")
figplt.set_xlabel("years")
figplt.set_ylabel("cohort ID")
Z=np.transpose(BA[pcid][pinst][1:150,0+pftadj:12+pftadj])
fig, figplt = plt.subplots(nrows=1,ncols=1, figsize=(20,40))
figplt.imshow(Z)
figplt.set_title("Basal area (m^2/ha)")
figplt.set_xlabel("years")
figplt.set_ylabel("cohort ID")
Z=np.transpose(NPLANT_CAN[pcid][pinst][1:150,0+pftadj:12+pftadj])
fig, figplt = plt.subplots(nrows=1,ncols=1, figsize=(20,40))
figplt.imshow(Z)
figplt.set_title("NPLANT CANOPY (N trees/ha)")
figplt.set_xlabel("years")
figplt.set_ylabel("cohort ID")
Z=np.transpose(DDBH_CANscpfadj[pcid][pinst][1:150*12,0+pftadj:12+pftadj])
fig, figplt = plt.subplots(nrows=1,ncols=1, figsize=(20,40))
figplt.imshow(Z)
figplt.set_title("DDBH of canopy layer (cm/yr/tree)")
figplt.set_xlabel("months")
figplt.set_ylabel("cohort ID")
Out[13]:
In [11]:
# plot recruitment
fig,figmet1 = plt.subplots(figsize=(10,3))
figmet1.plot(range(1,150-1), TBOTyrmean[1][1:150-1], c="red",label="temperature")
figmet1.set_ylabel("temperature")
figmet1.set_title("Mean annual TBOT and RAIN cz1")
figmet1b=figmet1.twinx()
figmet1b.plot(range(1,150-1), RAINyrmean[1][1:150-1], c="green",label="RAIN")
figmet1b.set_ylabel("RAIN (mm/s)")
figmet1.legend()
figmet1b.legend()
fig,figmet2 = plt.subplots(figsize=(10,3))
figmet2.plot(range(1,150-1), TBOTyrmean[4][1:150-1], c="red",label="temperature")
figmet2.set_ylabel("cz temperature")
figmet2.set_title("Mean annual TBOT and RAIN cz2")
figmet2b=figmet2.twinx()
figmet2b.plot(range(1,150-1), RAINyrmean[4][1:150-1], c="green",label="RAIN")
figmet2b.set_ylabel("RAIN (mm/s)")
figmet2.legend()
figmet2b.legend()
for pcid in cases:
for pinst in sets:
fig, (figplt1, figplt2, figplt3,figplt4) = plt.subplots(nrows=4,ncols=1, figsize=(10,8))
figplt1.plot(range(1,150*12), REC[pcid][pinst][1:150*12,0], linewidth=1,label="pine")
figplt1.plot(range(1,150*12), REC[pcid][pinst][1:150*12,1], linewidth=1,label="oak")
figplt1.set_title("recruitment:c"+str(pcid)+"ins"+str(pinst))
figplt1.set_xlabel("year")
figplt1.set_ylabel("Recruitment indiv/ha/yr")
figplt1.legend()
figplt2.plot(range(1,150-1), TBAyrmean[pcid][pinst][1:150-1], linewidth=1,c="blue",label="TBA")
figplt2.set_title("TBA and AWP:c"+str(pcid)+"ins"+str(pinst))
figplt2.set_xlabel("year")
figplt2.set_ylabel("TBA (m^2/ha)")
figplt2b=figplt2.twinx()
figplt2b.plot(range(1,150-1), AWPyrmean[pcid][pinst][1:150-1,13], c="red",label="AWP")
figplt2b.set_ylabel("AWP (MPa)")
figplt2.legend()
figplt2b.legend()
#figplt2c.legend()
figplt3.plot(range(1,150-1), TBAyrmean[pcid][pinst][1:150-1], linewidth=1,c="blue",label="TBA")
figplt3.set_title("TBA and TLAI:c"+str(pcid)+"ins"+str(pinst))
figplt3.set_xlabel("year")
figplt3.set_ylabel("TBA (m^2/ha)")
figplt3b=figplt3.twinx()
figplt3b.plot(range(1,150-1), TLAIyrmean[pcid][pinst][1:150-1], c="red",label="TLAI")
figplt3b.set_ylabel("LAI (m^2/m^2)")
figplt3.legend()
figplt3b.legend()
figplt4.plot(range(1,150-1), GPPyrmean[pcid][pinst][1:150-1], linewidth=1,c="blue",label="GPP")
figplt4.set_title("GPP and QVEGT:c"+str(pcid)+"ins"+str(pinst))
figplt4.set_xlabel("year")
figplt4.set_ylabel("GPP (fC/yr/ha)")
figplt4b=figplt4.twinx()
figplt4b.plot(range(1,150-1), QVEGTyrmean[pcid][pinst][1:150-1], c="red",label="QVEGT")
figplt4b.set_ylabel("QVEGT (gH2O/yr/ha)")
figplt4.legend()
figplt4b.legend()
fig, histpl = plt.subplots()
#histpl.hist(AWPyrmean[pcid][pinst][1:150-1,1+pftadj], bins=10)
In [18]:
# plot DDBHadj by month of canopy layer for each cohort bin
pcid=3
pftadj=13
#pinst=19
for pchid in range(0,12):
fig, figplt = plt.subplots(nrows=1,ncols=1, sharey=True, figsize=(8,2))
figplt.plot(range(1,150), DDBH_CANadjyrmean[pcid][pinst][1:150,pchid], linewidth=1,label='Canopy pine')
figplt.plot(range(1,150), DDBH_UNDadjyrmean[pcid][pinst][1:150,pchid], linewidth=1, label='Understory pine')
figplt.plot(range(1,150), DDBH_CANadjyrmean[pcid][pinst][1:150,pchid+pftadj], linewidth=1,label='Canopy oak')
figplt.plot(range(1,150), DDBH_UNDadjyrmean[pcid][pinst][1:150,pchid+pftadj], linewidth=1, label='Understory oak')
figplt.set_title("mean chort size (DBH cm):"+str(DBH[pchid]))
figplt.set_xlabel("year")
figplt.set_ylabel("DDBH canopy (cm/tree/yr)")
figplt.legend()
In [20]:
# plot DDBHadj by month of canopy layer for each cohort bin
pcid=3
pftadj=13
#pinst=19
for pchid in range(0,12):
fig, figplt = plt.subplots(nrows=1,ncols=1, sharey=True, figsize=(8,2))
figplt.plot(range(1,150), NPLANT_CANyrmean[pcid][pinst][1:150,pchid], linewidth=1,label='Canopy pine')
figplt.plot(range(1,150), NPLANTyrmean[pcid][pinst][1:150,pchid], linewidth=1, label='All pine')
figplt.plot(range(1,150), NPLANT_CANyrmean[pcid][pinst][1:150,pchid+pftadj], linewidth=1,label='Canopy oak')
figplt.plot(range(1,150), NPLANTyrmean[pcid][pinst][1:150,pchid+pftadj], linewidth=1, label='All oak')
figplt.set_title("Number plants (DBH cm):"+str(DBH[pchid]))
figplt.set_xlabel("year")
figplt.set_ylabel("Number (N tree/ha)")
figplt.legend()
In [16]:
# plot DDBHadj yrmean for each cohort bin
#pcid=5
#pinst=2
for pchid in range(0,12):
fig, figplt = plt.subplots(nrows=1,ncols=1, sharey=True, figsize=(8,2))
figplt.plot(range(1,150), 10000*DBA_CANadjyrmean[pcid][pinst][1:150,pchid+pftadj], linewidth=1,label='DBA')
# figplt.plot(range(1,150), 10000*MTba_CANadjyrmean[pcid][pinst][1:150,pchid], linewidth=1)
figplt.set_title("mean chort size of canopy layer (DBH cm):"+str(DBH[pchid]))
figplt.set_xlabel("year")
figplt.set_ylabel("cm^2/tree/yr")
figplt.legend()
print("growth and mortality as basal area")
In [17]:
# plot annual MT rate yrmean for each cohort bin
pcid=5
#pinst=2
for pchid in range(0,12):
fig, figplt = plt.subplots(nrows=1,ncols=1, sharey=True, figsize=(8,2))
# figplt.plot(range(1,150), DBA_UNDadjyrmean[pcid][pinst][1:150,pchid], linewidth=1)
figplt.plot(range(1,150*12), MT_CANscpfadj[pcid][pinst][1:150*12,pchid+pftadj], linewidth=1,label='Canopy')
figplt.plot(range(1,150*12), MTadj[pcid][pinst][1:150*12,pchid], linewidth=1,label='All')
figplt.set_title("mean chort size understory layer (DBH cm):"+str(DBH[pchid]))
figplt.set_xlabel("month")
figplt.set_ylabel("mortality rate (%/yr)")
figplt.legend()
In [18]:
# plot NPLANT for each cohort bin
#pcid=5
#pinst=20
for pchid in range(0,12):
fig, figplt = plt.subplots(nrows=1,ncols=1, sharey=True, figsize=(8,2))
figplt.plot(range(1,150*12), NPLANT[pcid][pinst][1:150*12,pchid+pftadj], linewidth=1)
figplt.plot(range(1,150*12), OTHER3scpf[pcid][pinst][1:150*12,pchid+pftadj], linewidth=1)
figplt.legend(["all", "canopy"])
figplt.set_title("mean chort size (DBH cm):"+str(DBH[pchid]))
figplt.set_xlabel("months")
figplt.set_ylabel("NPLANT (n/ha)")
In [19]:
# plot NPP_SCPF for each cohort bin
#pcid=5
#pinst=20
for pchid in range(0,12):
fig, figplt = plt.subplots(nrows=1,ncols=1, sharey=True, figsize=(8,2))
figplt.plot(range(1,150*12), NPPscpfadj[pcid][pinst][1:150*12,pchid+pftadj], linewidth=1)
figplt.set_title("mean chort size (DBH cm):"+str(DBH[pchid]))
figplt.set_xlabel("months")
figplt.set_ylabel("NPP per tree (kgC/year/tree)")
In [20]:
# plot AWP for each cohort bin
#pcid=5
#pinst=2
for pchid in range(0,12):
fig, figplt = plt.subplots(nrows=1,ncols=1, sharey=True, figsize=(8,2))
figplt.plot(range(1,150), AWPyrmean[pcid][pinst][1:150,pchid+pftadj], linewidth=1)
figplt.set_title("AWP (MPa):"+str(DBH[pchid]))
figplt.set_xlabel("year")
figplt.set_ylabel("AWP (MPa)")
fig, fig1plt = plt.subplots(nrows=1,ncols=1, sharey=True, figsize=(8,2))
for pchid in range(0,12):
fig1plt.plot(range(1,150), AWPyrmean[pcid][pinst][1:150,pchid+pftadj], linewidth=0.5, label=str(DBH[pchid]))
fig1plt.plot(range(1,150), AWPpftyrmean[pcid][pinst][1:150,pft], linewidth=1, label="all",c="red")
fig1plt.set_title("AWP (MPa)")
fig1plt.set_xlabel("year")
fig1plt.set_ylabel("AWP (MPa)")
fig1plt.legend()
Out[20]:
In [26]:
# define curve fitting functions for TBA vs AWP, MT and DBA vs AWP
# use linear function for TBA vs AWP
def fitTBA(awp,a,b):
tba=a+b*awp
return(tba)
# use power law function for MT and DBA
def fitMT(awp,inter,coe,exp):
mt=inter+coe*awp**exp
return(mt)
def fitDBA(awp,inter,coe,exp):
dba=inter+coe*awp**exp
return(dba)
In [29]:
# putup a master table for BA with rain
inset=np.array([5,11,29,35])
Thecid=4
Thestyr=1
Theendyr=149
Thechid=0
Thepft=PFT
# define arrays to contain the pars and covs from curve fitting for each instance
BAtb = pd.DataFrame({'rain':RAINyrmean[5][1:150-1]})
for insid in range(0,4):
cname="TBA"+str(inset[insid])
Cdata=TBApftyrmean[Thecid][inset[insid]][1:150-1,Thepft]
BAtb[cname]=Cdata.tolist()
cname="AWP"+str(inset[insid])
Cdata=AWPpftyrmean[Thecid][inset[insid]][1:150-1,Thepft]
BAtb[cname]=Cdata.tolist()
# df['blocks'] = blocks.tolist()
#BAtb["rn"]=[10, 20,30]
#df[df["A"] > 0]
RAINqtls=BAtb['rain'].quantile([0.05, 0.475, 0.525,0.95])
BAtbq10=BAtb[BAtb['rain']<=RAINqtls[0.05]]
BAtbq50=BAtb[(BAtb['rain']>=RAINqtls[0.475]) & (BAtb['rain']>=RAINqtls[0.525])]
BAtbq90=BAtb[BAtb['rain']>=RAINqtls[0.95]]
# fit curves to each rain range curve_fit(f, xdata, ydata)
#curve_fit(fitTBA, xdata, ydata)
# plot BApft vs. AWPpft
fig,TBAs=plt.subplots(nrows=2,ncols=2,figsize=(11,11))
insid=0
for i in range(0,2):
for j in range(0,2):
# plot DBA vs AWP
xname="AWP"+str(inset[insid])
yname="TBA"+str(inset[insid])
TBAs[i][j].scatter(AWPyrmean[Thecid][inset[insid]][Thestyr:Theendyr,Thechid+pftadj], TBApftyrmean[cid][inset[insid]][Thestyr:Theendyr,Thepft],c="gray") # ,marker=MK[1:54]
# fit linear function to each rain range, get estimated parameters
pars_q10, cov_q10 = curve_fit(fitTBA, BAtbq10[xname], BAtbq10[yname])
pars_q50, cov_q50 = curve_fit(fitTBA, BAtbq50[xname], BAtbq50[yname])
pars_q90, cov_q90 = curve_fit(fitTBA, BAtbq90[xname], BAtbq90[yname])
# plot raw data
TBAs[i][j].scatter(BAtbq10[xname],BAtbq10[yname],c="red")
TBAs[i][j].scatter(BAtbq50[xname],BAtbq50[yname],c="green")
TBAs[i][j].scatter(BAtbq90[xname],BAtbq90[yname],c="blue")
# plot fitted curve
Rangeq10=np.arange(min(BAtbq10[xname]), max(BAtbq10[xname]),0.1) # x range of dataset
Rangeq50=np.arange(min(BAtbq50[xname]), max(BAtbq50[xname]),0.1)
Rangeq90=np.arange(min(BAtbq90[xname]), max(BAtbq90[xname]),0.1)
#TBAs[i][j].plot(Rangeq10,fitTBA(Rangeq10,*pars_q10),c="red")
#TBAs[i][j].plot(Rangeq50,fitTBA(Rangeq50,*pars_q50),c="green")
#TBAs[i][j].plot(Rangeq90,fitTBA(Rangeq90,*pars_q90),c="blue")
TBAs[i][j].set_xlabel('AWP', labelpad = 5)
TBAs[i][j].set_ylabel('TBA (m^2/ha)', labelpad = 5)
TBAs[i][j].set_title("Mean annual TBA vs AWP, ins:"+str(inset[insid]))
# TBAs[i][j].set_ylim([6, 12])
insid=insid+1
# BAtb['rain'].hist
In [28]:
# plot DBA vs AWP and MTba vs AWP on the same plot
# Sept 9th working on this
inset=np.array([5,11,29,35])
#Thecid=5
Thestyr=1
Theendyr=148
Thechid=0
Thepft=2
# define arrays to contain par and covs from curve fitting
par_mtba = np.empty((len(path),nins+1,13*3,3))
cov_mtba = np.empty((len(path),nins+1,13*3,3,3))
par_dba = np.empty((len(path),nins+1,13*3,3))
cov_dba = np.empty((len(path),nins+1,13*3,3,3))
for cid in cases:
for sid in inset:
xvar = np.squeeze(AWPyrmean[cid][sid][Thestyr:Theendyr,Thechid+pftadj])
yvar_dba = np.squeeze(10000*DBA_CANadjyrmean[cid][sid][Thestyr:Theendyr,Thechid+pftadj])
yvar_mtba = np.squeeze(10000*MTba_CANadjyrmean[cid][sid][Thestyr:Theendyr,Thechid+pftadj])
#print("case:", cid, "set:", sid)
#print("xvar:",xvar)
#print("yvar_dba:",yvar_dba)
try:
pardba,covdba=curve_fit(fitDBA, xvar, yvar_dba)
parmt,covmt=curve_fit(fitMT, xvar, yvar_mtba)
par_dba[cid][sid][Thechid+pftadj]=pardba
#print("pardba",pardba)
cov_dba[cid][sid][Thechid+pftadj]=covdba
par_mtba[cid][sid][Thechid+pftadj]=parmt
cov_mtba[cid][sid][Thechid+pftadj]=covmt
except ValueError:
print("case:", cid, "set:", sid)
pass
print("Plot DBA and MT for canopy layer")
# plot DBA and MT of canopy layer vs TBOT
for Thecid in cases:
fig,BAs=plt.subplots(nrows=2,ncols=2,figsize=(11,11))
#plt.scatter(DifBAmeanLM1, DifBAmeanHM1)
insid=0
for i in range(0,2):
for j in range(0,2):
# plot DBA vs AWP
BAs[i][j].scatter(AWPyrmean[Thecid][inset[insid]][Thestyr:Theendyr,Thechid+pftadj], 10000*DBA_CANadjyrmean[cid][inset[insid]][Thestyr:Theendyr,Thechid+pftadj]) # ,marker=MK[1:54]
BAs[i][j].scatter(AWPyrmean[Thecid][inset[insid]][Thestyr:Theendyr,Thechid+pftadj], 10000*MTba_CANadjyrmean[cid][inset[insid]][Thestyr:Theendyr,Thechid+pftadj]) # ,marker=MK[1:54]
#Rangex=np.arange(min(AWPyrmean[Thecid][inset[insid]][Thestyr:Theendyr,Thechid+pftadj]), max(AWPyrmean[Thecid][inset[insid]][Thestyr:Theendyr,Thechid+pftadj]),0.1)
#BAs[i][j].plot(Rangex,fitDBA(Rangex,*par_dba[Thecid][inset[insid]][Thechid+pftadj][:]),c="blue")
#BAs[i][j].plot(Rangex,fitMT(Rangex,*par_dba[Thecid][inset[insid]][Thechid+pftadj][:]),c="red")
BAs[i][j].set_xlabel('AWP', labelpad = 5)
BAs[i][j].set_ylabel('cm^2/yr/tree', labelpad = 5)
BAs[i][j].set_title("Mean annual growth and mortality vs AWP, c"+str(Thecid)+"ins"+str(inset[insid]))
#BAs[i][j].set_ylim([0, 0.00003])
insid=insid+1
# TBOT
fig,BAs=plt.subplots(nrows=2,ncols=2,figsize=(11,11))
#plt.scatter(DifBAmeanLM1, DifBAmeanHM1)
insid=0
for i in range(0,2):
for j in range(0,2):
# plot DBA vs AWP
BAs[i][j].scatter(TBOTyrmean[Thecid][Thestyr:Theendyr], 10000*DBA_CANadjyrmean[cid][inset[insid]][Thestyr:Theendyr,Thechid+pftadj],label="DBA") # ,marker=MK[1:54]
BAs[i][j].scatter(TBOTyrmean[Thecid][Thestyr:Theendyr], 10000*MTba_CANadjyrmean[cid][inset[insid]][Thestyr:Theendyr,Thechid+pftadj],label="MTba") # ,marker=MK[1:54]
#Rangex=np.arange(min(AWPyrmean[Thecid][inset[insid]][Thestyr:Theendyr,Thechid+pftadj]), max(AWPyrmean[Thecid][inset[insid]][Thestyr:Theendyr,Thechid+pftadj]),0.1)
#BAs[i][j].plot(Rangex,fitDBA(Rangex,*par_dba[Thecid][inset[insid]][Thechid+pftadj][:]),c="blue")
#BAs[i][j].plot(Rangex,fitMT(Rangex,*par_dba[Thecid][inset[insid]][Thechid+pftadj][:]),c="red")
BAs[i][j].set_xlabel('TBOT K', labelpad = 5)
BAs[i][j].set_ylabel('cm^2/yr/tree', labelpad = 5)
BAs[i][j].set_title("Mean annual growth and mortality vs TBOT, c"+str(Thecid)+"ins"+str(inset[insid]))
BAs[i][j].legend()
insid=insid+1
fig,BAvsMT=plt.subplots(nrows=2,ncols=2,figsize=(11,11))
insid=0
for i in range(0,2):
for j in range(0,2):
# plot DBA vs MTba
BAvsMT[i][j].scatter(10000*MTba_CANadjyrmean[cid][inset[insid]][Thestyr:Theendyr,Thechid+pftadj],10000*DBA_CANadjyrmean[cid][inset[insid]][Thestyr:Theendyr,Thechid+pftadj]) # ,marker=MK[1:54]
# BAvsMT[i][j].scatter(AWPyrmean[Thecid][inset[insid]][Thestyr:Theendyr,Thechid-max(i,j)], MTbayrmean[cid][inset[insid]][Thestyr:Theendyr,Thechid-max(i,j)]) # ,marker=MK[1:54]
BAvsMT[i][j].set_xlabel('MTba (cm^2/year/tree)', labelpad = 5)
BAvsMT[i][j].set_ylabel('DBA (cm^2/yr/tree)', labelpad = 5)
BAvsMT[i][j].set_title("Growth vs mortality, ins:c"+str(Thecid)+"ins"+str(inset[insid]))
#BAvsMT[i][j].set_ylim([0, 0.00002])
#BAvsMT[i][j].set_xlim([0, 0.00002])
insid=insid+1
In [24]:
BAmean[2][20:38]
Out[24]:
BAmean[3][20:38]
DifBAmeanHM1[20:38]
In [40]:
# compute the difference in BAmean,std; DDBHmean,DDBHstd; MTmean, std between low-medium,
# high - medium co2 by per ppm change of c02
DifCO2_low_med=386-283
DifCO2_high_med=566-386
DifBAmeanLM1=(BAmean[2][:]-BAmean[1][:])/DifCO2_low_med # cz1
DifBAmeanHM1=(BAmean[3][:]-BAmean[2][:])/DifCO2_high_med # cz1
DifBAmeanLM2=(BAmean[5][:]-BAmean[4][:])/DifCO2_low_med # cz2
DifBAmeanHM2=(BAmean[6][:]-BAmean[5][:])/DifCO2_high_med # cz2
DifBAstdLM1=(BAstd[2][:]-BAstd[1][:])/DifCO2_low_med # cz1
DifBAstdHM1=(BAstd[3][:]-BAstd[2][:])/DifCO2_high_med # cz1
DifBAstdLM2=(BAstd[5][:]-BAstd[4][:])/DifCO2_low_med # cz2
DifBAstdHM2=(BAstd[6][:]-BAstd[5][:])/DifCO2_high_med # cz2
DifMTbameanLM1=(MTbamean[2][:]-MTbamean[1][:])/DifCO2_low_med # cz1
DifMTbameanHM1=(MTbamean[3][:]-MTbamean[2][:])/DifCO2_high_med # cz1
DifMTbameanLM2=(MTbamean[5][:]-MTbamean[4][:])/DifCO2_low_med # cz2
DifMTbameanHM2=(MTbamean[6][:]-MTbamean[5][:])/DifCO2_high_med # cz2
DifMTbastdLM1=(MTbastd[2][:]-MTbastd[1][:])/DifCO2_low_med # cz1
DifMTbastdHM1=(MTbastd[3][:]-MTbastd[2][:])/DifCO2_high_med # cz1
DifMTbastdLM2=(MTbastd[5][:]-MTbastd[4][:])/DifCO2_low_med # cz2
DifMTbastdHM2=(MTbastd[6][:]-MTbastd[5][:])/DifCO2_high_med # cz2
DifDBAmeanLM1=(DBAmean[2][:]-DBAmean[1][:])/DifCO2_low_med # cz1
DifDBAmeanHM1=(DBAmean[3][:]-DBAmean[2][:])/DifCO2_high_med # cz1
DifDBAmeanLM2=(DBAmean[5][:]-DBAmean[4][:])/DifCO2_low_med # cz2
DifDBAmeanHM2=(DBAmean[6][:]-DBAmean[5][:])/DifCO2_high_med # cz2
DifDBAstdLM1=(DBAstd[2][:]-DBAstd[1][:])/DifCO2_low_med # cz1
DifDBAstdHM1=(DBAstd[3][:]-DBAstd[2][:])/DifCO2_high_med # cz1
DifDBAstdLM2=(DBAstd[5][:]-DBAstd[4][:])/DifCO2_low_med # cz2
DifDBAstdHM2=(DBAstd[6][:]-DBAstd[5][:])/DifCO2_high_med # cz2
In [41]:
# PLOT
# compare dif l-m and h-m
# set the colors
clid=1
CR=[0.2,0.6,1]
CG=[0.4,0.8,1]
CB=[0.2,0.6,1]
Mark=[12,20,28] #["o","v","x","."]
CO2colors=np.zeros([55,3])
MK=np.full(55,0)
xx=np.zeros((2,2,2))
#defube instances to examine
#insSet=np.array([20,23,26,29,32,35, 38,41,44,47,50,53])
# uncomment the line below for all instances
insSet=np.arange(1,54)
#MK=['o', 'v', 'x', '.', 'o', 'v', 'x', '.', 'o', 'v', 'x', '.', 'o',
# 'v', 'x', '.', 'o', 'v', 'x', '.', 'o', 'v', 'x', '.', 'o', 'v',
# 'x', '.', 'o', 'v', 'x', '.', 'o', 'v', 'x', '.', 'o', 'v', 'x',
# '.', 'o', 'v', 'x', '.', 'o', 'v', 'x', '.', 'x', 'x', 'x', 'x',
# 'x']
for vid in range(0,3):
for pgsid in range(0,2):
for agsid in range(0,3):
for rid in range(0,3):
CO2colors[clid]=(CR[vid],CG[pgsid],CB[agsid])
MK[clid]=Mark[rid]
clid=clid+1
fig,BAs=plt.subplots(nrows=2,ncols=2,figsize=(11,11))
#plt.scatter(DifBAmeanLM1, DifBAmeanHM1)
BAs[0][0].scatter(DifBAmeanLM1[insSet], DifBAmeanHM1[insSet],c=CO2colors[insSet],s=MK[insSet]) # ,marker=MK[1:54]
BAs[0][0].set_xlabel('M-L', labelpad = 5)
BAs[0][0].set_ylabel('H-M', labelpad = 5)
BAs[0][0].set_title("CZ1 BAmean")
xx[0][0]=np.array([min(DifBAmeanLM1),max(DifBAmeanLM1)]) # adding 1-1 line
BAs[0][1].scatter(DifBAstdLM1[insSet], DifBAstdHM1[insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[0][1].set_xlabel('M-L', labelpad = 5)
BAs[0][1].set_ylabel('H-M', labelpad = 5)
BAs[0][1].set_title("CZ1 BAstd")
xx[0][1]=np.array([min(DifBAstdLM1),max(DifBAstdLM1)])
BAs[1][0].scatter(DifBAmeanLM2[insSet], DifBAmeanHM2[insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[1][0].set_xlabel('M-L', labelpad = 5)
BAs[1][0].set_ylabel('H-M', labelpad = 5)
BAs[1][0].set_title("CZ2 BAmean")
xx[1][0]=np.array([min(DifBAmeanLM2),max(DifBAmeanLM2)])
BAs[1][1].scatter(DifBAstdLM2[insSet], DifBAstdHM2[insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[1][1].set_xlabel('M-L', labelpad = 5)
BAs[1][1].set_ylabel('H-M', labelpad = 5)
BAs[1][1].set_title("CZ2 BAstd")
xx[1][1]=np.array([min(DifBAstdLM2),max(DifBAstdLM2)])
for i in range(0,2):
for j in range(0,2):
BAs[i][j].plot([xx[i][j][0],xx[i][j][1]],[xx[i][j][0],xx[i][j][1]],"r--")
plt.show
print("difference in BAmean and BAstd between different CO2 treatment standardized by the diffrence in CO2 concentration between the corresponding two treatments")
print("")
print("Impact of co2 on the mean MT and DBA")
CO2colors
print(xx[0][1][0]) plt.plot([xx[0][1][0],xx[0][1][1]],[xx[0][1][0],xx[0][1][1]],"r--")
In [42]:
# PLOT
# compare dif l-m and h-m difDBA vs. MTba
fig,BAs=plt.subplots(nrows=2,ncols=2,figsize=(11,11))
#plt.scatter(DifBAmeanLM1, DifBAmeanHM1)
BAs[0][0].scatter(DifMTbameanLM1[insSet], DifDBAmeanLM1[insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[0][0].set_xlabel('MT (m^2/ha/yr)', labelpad = 5)
BAs[0][0].set_ylabel('DBA (m^2/ha/yr)', labelpad = 5)
BAs[0][0].set_title("CZ1 M-L ")
BAs[0][1].scatter(DifMTbameanHM1[insSet], DifDBAmeanHM1[insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[0][1].set_xlabel('MT (m^2/ha/yr)', labelpad = 5)
BAs[0][1].set_ylabel('DBA (m^2/ha/yr)', labelpad = 5)
BAs[0][1].set_title("CZ1 H-M")
BAs[1][0].scatter(DifMTbameanLM2[insSet], DifDBAmeanLM2[insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[1][0].set_xlabel('MT (m^2/ha/yr)', labelpad = 5)
BAs[1][0].set_ylabel('DBA (m^2/ha/yr)', labelpad = 5)
BAs[1][0].set_title("CZ2 M-L ")
BAs[1][1].scatter(DifMTbameanHM2[insSet], DifDBAmeanHM2[insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[1][1].set_xlabel('MT (m^2/ha/yr)', labelpad = 5)
BAs[1][1].set_ylabel('DBA (m^2/ha/yr)', labelpad = 5)
BAs[1][1].set_title("CZ2 H-M")
for i in range(0,2):
for j in range(0,2):
BAs[i][j].set_xlim([-0.002, 0.007])
BAs[i][j].set_ylim([-0.002, 0.007])
BAs[i][j].plot([-0.005,0.005],[-0.005,0.005],"r--")
plt.show
print("difference in MTmean and DBAmean between different CO2 treatment standardized by the diffrence in CO2 concentration between the corresponding two treatments")
print("")
print("Impact of co2 on the mean of MT and DBA")
In [43]:
# PLOT
# compare dif l-m and h-m difDBA vs. MTba
fig,BAs=plt.subplots(nrows=2,ncols=2,figsize=(11,11))
#plt.scatter(DifBAmeanLM1, DifBAmeanHM1)
BAs[0][0].scatter(DifMTbastdLM1[insSet], DifDBAstdLM1[insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[0][0].set_xlabel('MT (m^2/ha/yr)', labelpad = 5)
BAs[0][0].set_ylabel('DBA (m^2/ha/yr)', labelpad = 5)
BAs[0][0].set_title("CZ1 M-L ")
BAs[0][1].scatter(DifMTbastdHM1[insSet], DifDBAstdHM1[insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[0][1].set_xlabel('MT (m^2/ha/yr)', labelpad = 5)
BAs[0][1].set_ylabel('DBA (m^2/ha/yr)', labelpad = 5)
BAs[0][1].set_title("CZ1 H-M")
BAs[1][0].scatter(DifMTbastdLM2[insSet], DifDBAstdLM2[insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[1][0].set_xlabel('MT (m^2/ha/yr)', labelpad = 5)
BAs[1][0].set_ylabel('DBA (m^2/ha/yr)', labelpad = 5)
BAs[1][0].set_title("CZ2 M-L ")
BAs[1][1].scatter(DifMTbastdHM2[insSet], DifDBAstdHM2[insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[1][1].set_xlabel('MT (m^2/ha/yr)', labelpad = 5)
BAs[1][1].set_ylabel('DBA (m^2/ha/yr)', labelpad = 5)
BAs[1][1].set_title("CZ2 H-M")
for i in range(0,2):
for j in range(0,2):
BAs[i][j].set_xlim([-0.003, 0.01])
BAs[i][j].set_ylim([-0.003, 0.01])
BAs[i][j].plot([-0.003,0.01],[-0.003,0.01],"r--")
plt.show
plt.show
print("difference in MTstd and DBAstd between different CO2 treatment standardized by the diffrence in CO2 concentration between the corresponding two treatments")
print("")
print("Impact of co2 level on the variation of growth and mortality rate")
In [44]:
# set up for plotting -
# select case and instance to compare/plot
caseSet = np.array([2])
insSet = np.array([1,3])
## dinguish pft by line type, instance by color
## assign line type, and line color
lc = np.array(['blue','brown','orange','pink','green','purple','olive','red','peachpuff']) # line color
ls = np.array(['-','--']) # line style
In [45]:
CO2colors[50]
Out[45]:
In [47]:
## Plot BAstd vs. BAmean for selected instances of intermediate CO2
## setup plots to make
caseSet=np.array([2,5])
insSet=np.array([23,32,41,50]) #20,23,26,29,32,35, 38,41,44,47,50,53 insid > 35 -> kmax=2, insid<=40 -> kmax=3
years=150
# define mark color and line type to distinguish instances
lc = np.array(['blue','brown','orange','pink','green','purple','olive','red','peachpuff','darkblue','black']) # line color
ls = np.array(['-','--']) # line style
# define mark colors
colors=np.array(['blue','brown','orange','pink','green','purple','olive','red','peachpuff','darkblue','black','brown']) # line color
TheMarkers=["$20$","$23$","$26$","$29$","$32$","$35$", "$38$","$41$","$44$","$47$","$50$","$53$"]
ilc=0
# plot BAstd vs. BAmean for 2 sites
fig,BAs=plt.subplots(nrows=2, ncols=2,figsize=(11,11))
BAs[0][0].scatter(BAmean[caseSet[0]][insSet], BAstd[caseSet[0]][insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[0][0].set_title("CZ1 (6.5m root)")
BAs[0][0].set_xlabel("BAmean (m^2/ha)")
BAs[0][0].set_ylabel("BAstd (m^2/ha)")
BAs[0][1].scatter(BAmean[caseSet[1]][insSet], BAstd[caseSet[1]][insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[0][1].set_title("CZ2 (6.5m root)")
BAs[0][1].set_xlabel("BAmean (m^2/ha)")
BAs[0][1].set_ylabel("BAstd (m^2/ha)")
# BAmean and BAstd cz2 vs. cz1
BAs[1][0].scatter(BAmean[caseSet[0]][insSet], BAmean[caseSet[1]][insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[1][0].set_title("BAmean (m^2/ha/yr, 6.5m root)")
BAs[1][0].set_xlabel("CZ1")
BAs[1][0].set_ylabel("CZ2)")
BAs[1][1].scatter(BAstd[caseSet[0]][insSet], BAstd[caseSet[1]][insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[1][1].set_title("BAstd (m^2/ha/yr. 6.5m root)")
BAs[1][1].set_xlabel("CZ1")
BAs[1][1].set_ylabel("CZ2")
for txt in insSet:
BAs[0][0].annotate(txt, (BAmean[caseSet[0]][txt], BAstd[caseSet[0]][txt]))
BAs[0][1].annotate(txt, (BAmean[caseSet[1]][txt], BAstd[caseSet[1]][txt]))
BAs[1][0].annotate(txt, (BAmean[caseSet[0]][txt], BAmean[caseSet[1]][txt]))
BAs[1][1].annotate(txt, (BAstd[caseSet[0]][txt], BAstd[caseSet[1]][txt]))
plt.show
##===================
# plot MTstd vs. MTmean for 2 sites
fig,MTs=plt.subplots(nrows=2, ncols=2,figsize=(11,11))
MTs[0][0].scatter(MTbamean[caseSet[0]][insSet], MTbastd[caseSet[0]][insSet],c=CO2colors[insSet],s=MK[insSet])
MTs[0][0].set_title("CZ1 (6.5m root)")
MTs[0][0].set_xlabel("MTmean (m^2/ha/yr)")
MTs[0][0].set_ylabel("MTstd (m^2/ha/yr)")
MTs[0][1].scatter(MTbamean[caseSet[1]][insSet], MTbastd[caseSet[1]][insSet],c=CO2colors[insSet],s=MK[insSet])
MTs[0][1].set_title("CZ2 (6.5m root)")
MTs[0][1].set_xlabel("MTmean (m^2/ha/yr)")
MTs[0][1].set_ylabel("MTstd (m^2/ha/yr)")
# MTmean and MTstd cz2 vs. cz1
MTs[1][0].scatter(MTbamean[caseSet[0]][insSet], MTbamean[caseSet[1]][insSet],c=CO2colors[insSet],s=MK[insSet])
MTs[1][0].set_title("MTmean (m^2/ha/yr, 6.5m root)")
MTs[1][0].set_xlabel("CZ1")
MTs[1][0].set_ylabel("CZ2)")
MTs[1][1].scatter(MTbastd[caseSet[0]][insSet], MTbastd[caseSet[1]][insSet],c=CO2colors[insSet],s=MK[insSet])
MTs[1][1].set_title("MTstd (m^2/ha/yr. 6.5m root)")
MTs[1][1].set_xlabel("CZ1")
MTs[1][1].set_ylabel("CZ2")
for txt in insSet:
MTs[0][0].annotate(txt, (MTbamean[caseSet[0]][txt], MTbastd[caseSet[0]][txt]))
MTs[0][1].annotate(txt, (MTbamean[caseSet[1]][txt], MTbastd[caseSet[1]][txt]))
MTs[1][0].annotate(txt, (MTbamean[caseSet[0]][txt], MTbamean[caseSet[1]][txt]))
MTs[1][1].annotate(txt, (MTbastd[caseSet[0]][txt], MTbastd[caseSet[1]][txt]))
plt.show
### plot DBAmean vs. DBAstd
fig,MTs=plt.subplots(nrows=2, ncols=2,figsize=(11,11))
MTs[0][0].scatter(DBAmean[caseSet[0]][insSet], DBAstd[caseSet[0]][insSet],c=CO2colors[insSet],s=MK[insSet],label=insSet)
MTs[0][0].set_title("CZ1 (6.5m root)")
MTs[0][0].set_xlabel("DBAmean (m^2/ha/yr)")
MTs[0][0].set_ylabel("DBAstd (m^2/ha/yr)")
MTs[0][1].scatter(DBAmean[5][insSet], DBAstd[caseSet[1]][insSet],c=CO2colors[insSet],s=MK[insSet],label=insSet)
MTs[0][1].set_title("CZ2 (6.5m root)")
MTs[0][1].set_xlabel("DBAmean (m^2/ha/yr)")
MTs[0][1].set_ylabel("DBAstd (m^2/ha/yr)")
# DBAmean and DBAstd cz2 vs. cz1
MTs[1][0].scatter(DBAmean[caseSet[0]][insSet], DBAmean[caseSet[1]][insSet],c=CO2colors[insSet],s=MK[insSet],label=insSet)
MTs[1][0].set_title("DBAmean (m^2/ha/yr, 6.5m root)")
MTs[1][0].set_xlabel("CZ1")
MTs[1][0].set_ylabel("CZ2)")
MTs[1][1].scatter(DBAstd[caseSet[0]][insSet], DBAstd[caseSet[1]][insSet],c=CO2colors[insSet],s=MK[insSet],label=insSet)
MTs[1][1].set_title("DBAstd (m^2/ha/yr. 6.5m root)")
MTs[1][1].set_xlabel("CZ1")
MTs[1][1].set_ylabel("CZ2")
for txt in insSet:
MTs[0][0].annotate(txt, (DBAmean[caseSet[0]][txt], DBAstd[caseSet[0]][txt]))
MTs[0][1].annotate(txt, (DBAmean[5][txt], DBAstd[caseSet[1]][txt]))
MTs[1][0].annotate(txt, (DBAmean[caseSet[0]][txt], DBAmean[caseSet[1]][txt]))
MTs[1][1].annotate(txt, (DBAstd[caseSet[0]][txt], DBAstd[caseSet[1]][txt]))
plt.show
Out[47]:
In [60]:
## Plot BAstd vs. BAmean for selected instances of intermediate CO2, shallow roots
## setup plots to make
caseSet=np.array([2,5])
insSet=np.array([23,32,41,50]) #21,24,27,30,33,36, 39,42,45,48,51,54 insid > 35 -> kmax=2, insid<=40 -> kmax=3
years=150
# define mark color and line type to distinguish instances
lc = np.array(['blue','brown','orange','pink','green','purple','olive','red','peachpuff','darkblue','black']) # line color
ls = np.array(['-','--']) # line style
# define mark colors
colors=np.array(['blue','brown','orange','pink','green','purple','olive','red','peachpuff','darkblue','black','brown']) # line color
TheMarkers=["$21$","$24$","$27$","$30$","$33$","$36$", "$39$","$42$","$45$","$48$","$51$","$54$"]
ilc=0
# plot BAstd vs. BAmean for 2 sites
fig,BAs=plt.subplots(nrows=2, ncols=2,figsize=(11,11))
BAs[0][0].scatter(BAmean[2][insSet], BAstd[2][insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[0][0].set_title("CZ1 (1.5m root)")
BAs[0][0].set_xlabel("BAmean (m^2/ha)")
BAs[0][0].set_ylabel("BAstd (m^2/ha)")
BAs[0][1].scatter(BAmean[5][insSet], BAstd[5][insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[0][1].set_title("CZ2 (1.5m root)")
BAs[0][1].set_xlabel("BAmean (m^2/ha)")
BAs[0][1].set_ylabel("BAstd (m^2/ha)")
# BAmean and BAstd cz2 vs. cz1
BAs[1][0].scatter(BAmean[2][insSet], BAmean[5][insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[1][0].set_title("BAmean (m^2/ha, 1.5m root)")
BAs[1][0].set_xlabel("CZ1")
BAs[1][0].set_ylabel("CZ2)")
BAs[1][1].scatter(BAstd[2][insSet], BAstd[5][insSet],c=CO2colors[insSet],s=MK[insSet])
BAs[1][1].set_title("BAstd (m^2/ha. 1.5m root)")
BAs[1][1].set_xlabel("CZ1")
BAs[1][1].set_ylabel("CZ2")
for txt in insSet:
BAs[0][0].annotate(txt, (BAmean[2][txt], BAstd[2][txt]))
BAs[0][1].annotate(txt, (BAmean[5][txt], BAstd[5][txt]))
BAs[1][0].annotate(txt, (BAmean[2][txt], BAmean[5][txt]))
BAs[1][1].annotate(txt, (BAstd[2][txt], BAstd[5][txt]))
plt.show
##===================
# plot MTstd vs. MTmean for 2 sites
fig,MTs=plt.subplots(nrows=2, ncols=2,figsize=(11,11))
MTs[0][0].scatter(MTbamean[2][insSet], MTbastd[2][insSet],c=CO2colors[insSet],s=MK[insSet])
MTs[0][0].set_title("CZ1 (6.5m root)")
MTs[0][0].set_xlabel("MTmean (m^2/ha/yr)")
MTs[0][0].set_ylabel("MTstd (m^2/ha/yr)")
MTs[0][1].scatter(MTbamean[5][insSet], MTbastd[5][insSet],c=CO2colors[insSet],s=MK[insSet])
MTs[0][1].set_title("CZ2 (6.5m root)")
MTs[0][1].set_xlabel("MTmean (m^2/ha/yr)")
MTs[0][1].set_ylabel("MTstd (m^2/ha/yr)")
# MTmean and MTstd cz2 vs. cz1
MTs[1][0].scatter(MTbamean[2][insSet], MTbamean[5][insSet],c=CO2colors[insSet],s=MK[insSet])
MTs[1][0].set_title("MTmean (m^2/ha/yr, 6.5m root)")
MTs[1][0].set_xlabel("CZ1")
MTs[1][0].set_ylabel("CZ2)")
MTs[1][1].scatter(MTbastd[2][insSet], MTbastd[5][insSet],c=CO2colors[insSet],s=MK[insSet])
MTs[1][1].set_title("MTstd (m^2/ha/yr. 6.5m root)")
MTs[1][1].set_xlabel("CZ1")
MTs[1][1].set_ylabel("CZ2")
for txt in insSet:
MTs[0][0].annotate(txt, (MTbamean[2][txt], MTbastd[2][txt]))
MTs[0][1].annotate(txt, (MTbamean[5][txt], MTbastd[5][txt]))
MTs[1][0].annotate(txt, (MTbamean[2][txt], MTbamean[5][txt]))
MTs[1][1].annotate(txt, (MTbastd[2][txt], MTbastd[5][txt]))
plt.show
### plot DBAmean vs. DBAstd
fig,MTs=plt.subplots(nrows=2, ncols=2,figsize=(11,11))
MTs[0][0].scatter(DBAmean[2][insSet], DBAstd[2][insSet],c=CO2colors[insSet],s=MK[insSet],label=insSet)
MTs[0][0].set_title("CZ1 (1.5m root)")
MTs[0][0].set_xlabel("DBAmean (m^2/ha/yr)")
MTs[0][0].set_ylabel("DBAstd (m^2/ha/yr)")
MTs[0][1].scatter(DBAmean[5][insSet], DBAstd[5][insSet],c=CO2colors[insSet],s=MK[insSet],label=insSet)
MTs[0][1].set_title("CZ2 (1.5m root)")
MTs[0][1].set_xlabel("DBAmean (m^2/ha/yr)")
MTs[0][1].set_ylabel("DBAstd (m^2/ha/yr)")
# DBAmean and DBAstd cz2 vs. cz1
MTs[1][0].scatter(DBAmean[2][insSet], DBAmean[5][insSet],c=CO2colors[insSet],s=MK[insSet],label=insSet)
MTs[1][0].set_title("DBAmean (m^2/ha/yr, 1.5m root)")
MTs[1][0].set_xlabel("CZ1")
MTs[1][0].set_ylabel("CZ2)")
MTs[1][1].scatter(DBAstd[2][insSet], DBAstd[5][insSet],c=CO2colors[insSet],s=MK[insSet],label=insSet)
MTs[1][1].set_title("DBAstd (m^2/ha/yr. 1.5m root)")
MTs[1][1].set_xlabel("CZ1")
MTs[1][1].set_ylabel("CZ2")
for txt in insSet:
MTs[0][0].annotate(txt, (DBAmean[2][txt], DBAstd[2][txt]))
MTs[0][1].annotate(txt, (DBAmean[5][txt], DBAstd[5][txt]))
MTs[1][0].annotate(txt, (DBAmean[2][txt], DBAmean[5][txt]))
MTs[1][1].annotate(txt, (DBAstd[2][txt], DBAstd[5][txt]))
plt.show
print("BAmean and BAstd for shallow roots, intermediate co2")
In [61]:
MTbamean[2][insSet]
Out[61]:
In [62]:
DBAmean[2][insSet]
Out[62]:
In [48]:
## Plot DDBH, REC, vs MTba for selected instances of intermediate CO2
## setup plots to make
caseSet=np.array([2,5])
insSet=np.array([23,32,41,50]) #20,23,26,29,32,35, 38,41,44,47,50,53 insid > 35 -> kmax=2, insid<=40 -> kmax=3
# define mark color and line type to distinguish instances
lc = np.array(['blue','brown','orange','pink','green','purple','olive','red','peachpuff']) # line color
ls = np.array(['-','--']) # line style
ilc=0
# plot total DDBH vs. MTba
fig,(Bt1,Bt2)=plt.subplots(nrows=1,ncols=2,figsize=(11,5))
Bt1.scatter(MTbamean[2][insSet], DBAmean[2][insSet],c=CO2colors[insSet],s=MK[insSet])
Bt2.scatter(MTbamean[5][insSet], DBAmean[5][insSet],c=CO2colors[insSet],s=MK[insSet])
Bt1.set_title("CZ1 TBA (6.5m root)")
Bt1.set_xlabel("Total Mortality as BA")
Bt1.set_ylabel("DBA (m^2/ha/yr)")
#Bt1.legend(insSet, title="Instance",loc='upper left')
Bt2.set_title("CZ2 TBA (6.5m root)")
Bt2.set_xlabel("Total Mortality as BA (m^2/ha/yr)")
Bt2.set_ylabel("DBA (m^2/ha/yr)")
#Bt2.legend(insSet, title="Instance",loc='upper left')
for txt in insSet:
Bt1.annotate(txt, (MTbamean[2][txt], DBAmean[2][txt]))
Bt2.annotate(txt, (MTbamean[5][txt], DBAmean[5][txt]))
plt.show
# plot recruitment vs. MTba
# plot the total DDBH + recruitment vs. MTba
fig,(Bt1,Bt2)=plt.subplots(nrows=1,ncols=2,figsize=(11,5))
Bt1.scatter((MTbamean[2][insSet]+RECbamean[2][insSet]), DBAmean[2][insSet],c=CO2colors[insSet],s=MK[insSet])
Bt2.scatter((MTbamean[5][insSet]+RECbamean[5][insSet]), DBAmean[5][insSet],c=CO2colors[insSet],s=MK[insSet])
Bt1.set_title("CZ1 TBA (6.5m root)")
Bt1.set_xlabel("Total Mortality as BA")
Bt1.set_ylabel("DBA+Rec as BA (m^2/ha/yr)")
#Bt1.legend(insSet, title="Instance",loc='upper left')
Bt2.set_title("CZ2 TBA (6.5m root)")
Bt2.set_xlabel("Total Mortality as BA (m^2/ha/yr)")
Bt2.set_ylabel("DBA+Rec as BA (m^2/ha/yr)")
#Bt2.legend(insSet, title="Instance",loc='upper left')
plt.show
print("Total BA of Pine with 6.5m roots")
In [49]:
## Plot DDBH, REC, vs MTba for selected instances of intermediate CO2
## setup plots to make
caseSet=np.array([2,5])
insSet=np.array([23,32,41,50]) #21,24,27,30,33,36, 39,42,45,48,51,54 insid > 35 -> kmax=2, insid<=40 -> kmax=3
# define mark color and line type to distinguish instances
lc = np.array(['blue','brown','orange','pink','green','purple','olive','red','peachpuff']) # line color
ls = np.array(['-','--']) # line style
ilc=0
# plot total DDBH vs. MTba
fig,(Bt1,Bt2)=plt.subplots(nrows=1,ncols=2,figsize=(11,5))
Bt1.scatter(MTbamean[2][insSet], DBAmean[2][insSet],c=CO2colors[insSet],s=MK[insSet])
Bt2.scatter(MTbamean[5][insSet], DBAmean[5][insSet],c=CO2colors[insSet],s=MK[insSet])
Bt1.set_title("CZ1 TBA (1.5m root)")
Bt1.set_xlabel("Total Mortality as BA")
Bt1.set_ylabel("DBA (m^2/ha/yr)")
#Bt1.legend(insSet, title="Instance",loc='upper left')
Bt2.set_title("CZ2 TBA (1.5m root)")
Bt2.set_xlabel("Total Mortality as BA (m^2/ha/yr)")
Bt2.set_ylabel("DBA (m^2/ha/yr)")
#Bt2.legend(insSet, title="Instance",loc='upper left')
for txt in insSet:
Bt1.annotate(txt, (MTbamean[2][txt], DBAmean[2][txt]))
Bt2.annotate(txt, (MTbamean[5][txt], DBAmean[5][txt]))
plt.show
# plot recruitment vs. MTba
# plot the total DDBH + recruitment vs. MTba
fig,(Bt1,Bt2)=plt.subplots(nrows=1,ncols=2,figsize=(11,5))
Bt1.scatter((MTbamean[2][insSet]+RECbamean[2][insSet]), DBAmean[2][insSet],c=CO2colors[insSet],s=MK[insSet])
Bt2.scatter((MTbamean[5][insSet]+RECbamean[5][insSet]), DBAmean[5][insSet],c=CO2colors[insSet],s=MK[insSet])
Bt1.set_title("CZ1 TBA (1.5m root)")
Bt1.set_xlabel("Total Mortality as BA")
Bt1.set_ylabel("DBA+Rec as BA (m^2/ha/yr)")
#Bt1.legend(insSet, title="Instance",loc='upper left')
Bt2.set_title("CZ2 TBA (1.5m root)")
Bt2.set_xlabel("Total Mortality as BA (m^2/ha/yr)")
Bt2.set_ylabel("DBA+Rec as BA (m^2/ha/yr)")
#Bt2.legend(insSet, title="Instance",loc='upper left')
plt.show
print(" DBA vs MT of 1.5m roots instances")
In [56]:
## exploraty plots
## plot time series BA in different ways: A - compare BA between instance for the same case for 130 years
## plot different instance but same CO2 treatment in same subplot for each site
## setup plots to make
caseSet=np.array([2,5])
insSet=np.array([23,32,41,50])
#insSet=np.arange(1,54) 20,23,26,29,32,35, 38,41,44,47,50,53
#
# insid > 35 -> kmax=2, insid<=40 -> kmax=3
years=150
# define line color and line type to distinguish instances
lc = CO2colors #np.array(['blue','brown','orange','pink','green','purple','olive','red','peachpuff']) # line color
ls = np.array(['-','--']) # line style
ilc=0
fig,(Bt1,Bt2)=plt.subplots(nrows=2,ncols=1,figsize=(13,10))
for inst in insSet:
for cid in caseSet:
if inst <= 35:
ils = 0
else:
ils = 1
if cid <= 3:
Bt1.plot(Thetime[1:12*years], BAtotal[cid][inst][1:12*years], linestyle=ls[ils], color=lc[inst], linewidth=1)
else:
Bt2.plot(Thetime[1:12*years], BAtotal[cid][inst][1:12*years], linestyle=ls[ils], color=lc[inst], linewidth=1)
ilc=ilc+1
if ilc>=len(insSet)/2:
ilc=0
Bt1.set_title("CZ1 TBA (6.5m root)")
Bt1.set_xlabel("day since 1980-01-01")
Bt1.set_ylabel("BAtotal (m^2/ha)")
Bt1.legend(insSet, title="Instance",loc='upper left')
Bt2.set_title("CZ2 TBA (6.5m root)")
Bt2.set_xlabel("day since 1980-01-01")
Bt2.set_ylabel("BAtotal (m^2/ha)")
Bt2.legend(insSet, title="Instance",loc='upper left')
plt.show
print("Total BA of Pine with 6.5m roots, kmax=3,P50x=, differ in stoma")
In [66]:
## exploraty plots
## plot time series DDBHba in different ways: A - compare BA between instance for the same case for 130 years
## plot different instance but same CO2 treatment in same subplot for each site
## setup plots to make
DDBHtoMT=DDBHtotal-MTbatotal
caseSet=np.array([2,5])
insSet=np.array([20,23,26,29,32,35, 38,41,44,47,50,53]) # insid > 35 -> kmax=2, insid<=40 -> kmax=3
years=150
# define line color and line type to distinguish instances
lc = np.array(['blue','brown','orange','pink','green','purple','olive','red','peachpuff']) # line color
ls = np.array(['-','--']) # line style
ilc=0
#fig,(Bt1,Bt2)=plt.subplots(nrows=2,ncols=1,figsize=(13,13))
for inst in insSet:
fig,(Bt1,Bt2)=plt.subplots(nrows=2,ncols=1,figsize=(13,5))
for cid in caseSet:
if inst <= 35:
ils = 0
else:
ils = 1
if cid <= 3:
Bt1.plot(Thetime[1:12*years], DDBHtoMT[cid][inst][1:12*years], linestyle=ls[ils], color=lc[ilc], linewidth=1)
else:
Bt2.plot(Thetime[1:12*years], DDBHtoMT[cid][inst][1:12*years], linestyle=ls[ils], color=lc[ilc], linewidth=1)
ilc=ilc+1
if ilc>=len(insSet)/2:
ilc=0
Bt1.set_title("CZ1 Total DDBH-MTba (6.5m root), inst: "+str(inst),y=0.8)
Bt1.set_xlabel("day since 1980-01-01")
Bt1.set_ylabel("Total DDBH-MTba (m^2/ha/YR)")
# Bt1.legend(inst, title="Instance",loc='upper left')
Bt2.set_title("CZ2 Total DDBH-MTba (6.5m root), inst: "+str(inst), y=0.8)
Bt2.set_xlabel("day since 1980-01-01")
Bt2.set_ylabel("Total DDBH-MTba (m^2/ha/yr)")
# Bt2.legend(inst, title="Instance",loc='upper left')
plt.show
print("Total BA of Pine with 6.5m roots, kmax=3,P50x=, differ in stoma")
In [ ]:
## exploraty plots
## plot time series Recruitment between instance for the same case for 130 years
## plot different instance but same CO2 treatment in same subplot for each site
## setup plots to make
caseSet=np.array([2,5])
insSet=np.array([20,23,26,29,32,35,38,41,44,47,50,53]) # insid > 35 -> kmax=2, insid<=40 -> kmax=3
years=150
# define line color and line type to distinguish instances
lc = np.array(['blue','brown','orange','pink','green','purple','olive','red','peachpuff']) # line color
ls = np.array(['-','--']) # line style
ilc=0
fig,(Bt1,Bt2)=plt.subplots(nrows=2,ncols=1,figsize=(13,13))
for inst in insSet:
for cid in caseSet:
if inst <= 35:
ils = 0
else:
ils = 1
if cid <= 3:
Bt1.plot(Thetime[1:12*years], REC[cid][inst][1:12*years], linestyle=ls[ils], color=lc[ilc], linewidth=1)
else:
Bt2.plot(Thetime[1:12*years], REC[cid][inst][1:12*years], linestyle=ls[ils], color=lc[ilc], linewidth=1)
ilc=ilc+1
if ilc>=len(insSet)/2:
ilc=0
Bt1.set_title("CZ1 TBA (6.5m root)")
Bt1.set_xlabel("day since 1980-01-01")
Bt1.set_ylabel("Recruitment (n/ha/yr)")
Bt1.legend(insSet, title="Instance",loc='upper left')
Bt2.set_title("CZ2 TBA (6.5m root)")
Bt2.set_xlabel("day since 1980-01-01")
Bt2.set_ylabel("Recruitment (n/ha/yr)")
Bt2.legend(insSet, title="Instance",loc='upper left')
plt.show
print("Total BA of Pine with 6.5m roots, kmax=3,P50x=, differ in stoma")
In [ ]:
## exploraty plots
## plot time series BA in different ways: A - compare BA between instance for the same case for 130 years
## plot different CO2 treatment of same instances of CZ1 in same subplot for up to three instances
## setup plots to make
caseSet1=np.array([1,2,3]) # rerun case 3 - high co2
insSet1=np.array([20,23,26,29,32,35,38,41,44,47,50,53])
# shallow roots[]
years=150
# define line color and line type to distinguish instances
lc = np.array(['blue','brown','orange','pink','green','purple','olive','red','peachpuff']) # line color
ls = np.array(['-','--']) # line style
ils=0
fig,Bt1=plt.subplots(nrows=len(insSet1),ncols=1,figsize=(13,36))
for instid in range (0,len(insSet1)):
inst=insSet1[instid]
for cid in caseSet1:
print(cid,inst)
Bt1[instid].plot(Thetime[1:12*years], BAtotal[cid][inst][1:12*years], linestyle=ls[ils], color=lc[cid-1], linewidth=1)
Bt1[instid].set_title("CZ1 TBA ins:"+str(insSet1[instid]),y=0.1)
Bt1[instid].set_xlabel("day since 1980-01-01")
Bt1[instid].set_ylabel("BAtotal (m^2/ha)")
Bt1[instid].legend(caseSet1, title="CO2 level",loc='upper left')
plt.show
print("notes: the high c02 treatment drop to zero at arund 38000 days because the run stopped")
print("at that time due to time off, not because all the trees are dies suddently")
print("Also, case6 high co2 treatment were not run properly, running it again")
In [53]:
## exploraty plots
## plot time series BA in different ways: A - compare BA between instance for the same case for 130 years
## plot different CO2 treatment of same instances of CZ2 in same subplot for up to three instances
## setup plots to make
caseSet2=np.array([4,5,6]) # need to run high co2 case6 again
insSet2=np.array([20,23,26,29,32,35,38,41,44,47,50,53])
years=150
# define line color and line type to distinguish instances
lc = np.array(['blue','brown','orange','pink','green','purple','olive','red','peachpuff']) # line color
ls = np.array(['-','--']) # line style
ils=0
fig,Btx=plt.subplots(nrows=len(insSet2),ncols=1,figsize=(13,36))
for instid in range (0,len(insSet2)):
inst=insSet2[instid]
for cid in caseSet2:
print(cid,inst)
Btx[instid].plot(Thetime[1:12*years], BAtotal[cid][inst][1:12*years], linestyle=ls[ils], color=lc[cid-1], linewidth=1)
Btx[instid].set_title("CZ2 TBA ins:"+str(insSet2[instid]),y=0.1)
Btx[instid].set_xlabel("day since 1980-01-01")
Btx[instid].set_ylabel("BAtotal (m^2/ha)")
Btx[instid].legend(caseSet2, title="CO2 level",loc='upper left')
plt.show
print("notes: ")
In [ ]:
## plot GPP and QVEGT
## setup plots to make
# r1 long term GPP, r2 2yr GPP with flux tower, r3, long term qvegt, r4, 2yr LH with flux tower
sToyr=365*24*60*60
fig,(PGPP, PGPPflx,PQVEGT,PLHflx)=plt.subplots(nrows=4,ncols=1,figsize=(8,10))
fig.tight_layout()
for cid in caseSet:
for inst in insSet:
if inst <=8:
ils = 0
ilc = inst -1
else:
ils = 1
ilc = inst - 9
PGPP.plot(Thetime[1:len(Thetime)], GPP[inst][1:(len(Thetime))]*sToyr, linestyle=ls[ils], color=lc[ilc], linewidth=1)
PGPP.set_xlabel("day since 1980-01-01")
PGPP.set_ylabel("GPP (gC/m^2/yr)")
PGPP.legend(cases, title="Instance",loc='right')
# plot 2011 + 33 to 2012 + 33
st = (2066 - 1980) * 12-1
ent = st + 24
for inst in cases:
if inst <=8:
ils = 0
ilc = inst -1
else:
ils = 1
ilc = inst - 9
PGPPflx.plot(GPP[inst][st:ent]*sToyr, linestyle=ls[ils], color=lc[ilc],linewidth=1)
#PGPPflx.plot(cz1GEE*sToyr, linestyle='--', linewidth=2, color='blue')
PGPPflx.set_xlabel("month since 2011-01")
PGPPflx.set_ylabel("GPP (gC/m^2/yr)")
PGPPflx.legend(np.append(cases,'observed'), title="Instance", loc='right')
for inst in insSet:
if inst <=8:
ils = 0
ilc = inst -1
else:
ils = 1
ilc = inst - 9
PQVEGT.plot(Thetime[1:len(Thetime)], QVEGT[inst][1:(len(Thetime))], linestyle=ls[ils], color=lc[ilc],linewidth=1)
PQVEGT.set_xlabel("day since 1980-01-01")
PQVEGT.set_ylabel("QVEGT (mm/s)")
PQVEGT.legend(cases, title="Instance", loc='right')
for inst in insSet:
if inst <=8:
ils = 0
ilc = inst - 1
else:
ils = 1
ilc = inst - 9
PLHflx.plot(LH[inst][st:ent], linestyle=ls[ils], color=lc[ilc],linewidth=1)
#PLHflx.plot(cz1LH*sToyr, linestyle='--', linewidth=2, color='blue')
PLHflx.set_xlabel("month since 2011-01")
PLHflx.set_ylabel("Latent heat flux (W/m^2)")
PLHflx.ticklabel_format(useOffset=False)
PLHflx.legend(np.append(cases,'observed'), title="Instance", loc='right')
#+ 'observed'
print("notes: hightly over estimted WUE, this is likely due to higher")
print("oak proportion, oak has safer xylem and higher WUE")
print("solid line - fire off, dash line - fire on")
#plt.show
In [9]:
# plot MR to AR ratio
GRtoMR=GR/MR
GRtoMR.shape
Out[9]:
In [ ]:
# plot GR to AR ratio of size 55 cm DBH cohort of each pft
fig1,(Ppine,Poak) = plt.subplots(nrows=2,ncols=1,figsize=(10,8))
fig1.tight_layout()
for ins in range(1,nins+1):
Ppine.plot(Thetime[(len(Thetime)-240):len(Thetime)], GRtoMR[ins,(len(Thetime)-240):(len(Thetime)),8],linestyle='-', color=lc[ins-1])
Ppine.set_title("GR to MR ratio of :dbh 55cm")
Ppine.set_ylabel("GR/MR pine")
Ppine.set_xlabel("day since 1980-01-01")
Ppine.legend(range(1,nins+1), title="Instance",loc='upper left')
for ins in range(1,nins+1):
Poak.plot(Thetime[(len(Thetime)-240):len(Thetime)], GRtoMR[ins,(len(Thetime)-240):(len(Thetime)),21], linestyle='-', color=lc[ins-1])
#Poak.set_title("GR to MR ratio of Pine:dbh 55cm")
Poak.set_ylabel("GR/MR oak")
Poak.set_xlabel("day since 1980-01-01")
Poak.legend(range(1,nins+1), title="Instance",loc='upper left')
In [30]:
# plot BA of cohorts
fig1,pBAch = plt.subplots(nrows=8,figsize=(10,20))
fig1.tight_layout()
levels = np.arange(0.,1.1, 0.1)
cmap = plt.get_cmap('Blues')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)
for ins in range(1,nins+1):
im = pBAch[ins].pcolormesh(Thetime[1:len(Thetime)], Thesizebin, BA[ins,0:12], cmap=cmap, norm=norm)
pBAch[ins].set_title("GR to MR ratio of :dbh 55cm")
pBAch[ins].set_ylabel("Cohort ")
pBAch[ins].set_xlabel("day since 1980-01-01")
## set up the first plot: the fractional area of patches of a given age range
#levels = np.arange(0.,1.1, 0.1)
#cmap = plt.get_cmap('Blues')
#norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)
#im = f1ax0.pcolormesh(time, patch_age_bins, PATCH_AREA_BY_AGE[:,:,0].transpose(), cmap=cmap, norm=norm)
#fig1.colorbar(im, ax=f1ax0)
#f1ax0.set_title(r'Patch Area by Age (m$^2$ patch m$^{-2}$ site)')
#f1ax0.set_xlabel('Time (yr)')
#f1ax0.set_ylabel('Patch Age (yr)')
In [ ]: